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Appealing  to  current  and  past  work  on  the  routing  problem  in 
data  - communication  networks,  we  motivate  the  need  for  algor- 
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algorithms  that  process  the  record  of  arrivals  and  departures 
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ulated per  unit  time.  Through  simulation  and  analysis  we  show 
that  all  three  algorithms  are  asymptotically  unbiased  and 
efficient  for  M/D/1  queues.  By  simulation  of  other  queues  we 
investigate  the  relative  robustness  of  the  three  procedures. 
Finally,  through  examination  of  the  storage  and  computational 
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S1-;CT10N  1 
INTRODUCTION 


1.1  The  MessaRe-Routinp;  Problem  in  DaLa-Comniunlcati  Networks 

There  are  several  major  analytical  prohlcmis  in  the  design 
of  a modern  data  communication  network.  Given  a set  of  nodes, 
different  topological  configurations  may  be  considered.  Having 
specified  the  manner  in  which  the  nodes  are  connected,  there 
remains  the  question  of  how  to  assign  capacities  to  the  communi- 
cation links.  For  each  problem,  different  constraints  and 
optimization  criteria  are  appropriate.  Once  we  have  resolved 
the  issue  of  the  structure  of  the  network,  the  central  problem 
that  remains  is  how  to  route  a given  message  from  a source  to  a 
destination  node.  We  are  interested  in  s tore-and- forward  computer 
networks  where  messages,  or  segments  of  messages  called  packets, 
travel  from  a source  to  a destination  node,  waiting  in  queues 
for  retransmission  at  each  intermediate  node.  One  way  to  specify 
a routing  policy  is  by  providing  a routing  table  for  each  node  i 
listing  what  fraction  of  the  traffic  destined  for  node  j is  to 
be  sent  on  each  of  the  outgoitig  links  from  node  i.  Routing 
strategies  vary  in  character  from  purely  static,  the  routing 
fractions  being  fixed  in  time ^ and  determined  on  the  basis  of 
average  arrival  statistics,  to  the  completely  dynamic  case  where 
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Llic  frac- Lions  vary  conLinuousJy  in  time  according  Lo  Lhc  "staLc" 
of  Ulc  network.  The  philosophy  for  implementing  any  strategy 
varies  between  centralized  and  decentralized  extremes.  In  the 
centralized  case  a special  node  in  the  network  receives  informa- 
tion and  does  all  tlie  routing  computations,  communicating  changes 
in  the  routing  fractions  to  all  other  nodes.  In  a decentralized 
scheme,  each  nod  ■ computes  its  own  routing  table  on  the  basis  of 
"locally"  availa  le  information. 

An  intermed;  ary  between  the  strictly  static  c.nd  dynamic 
routing  schemes  as  termed  quasi-static  routing.  Here  we  overcome 
the  static  procedures  insensitivity  to  gradual  traffic  changes 
and  the  failures  of  links  and  nodes  by  up-dating  the  routi  g 
strategy  periodically,  or  when  a special  need  arises.  Ir.  the 
dynamic  routing  case,  messages  that  have  been  segmented  into 
packets  may  arrive  out  of  order  at  the  destination  node,  neces- 
sitating a "reassembly"  operation.  In  the  quasi-static  pro- 
cedure, most  messages  would  be  delivered  in  order  since  the  time 
intervals  between  routing  charges  will  be  relatively  long.  Hence, 
quasi-static  routing  procedures  suggest  a sensible  mid-point 
between  the  static  and  dynamic  extremes. 


1 . 2 of  llu-  Tliesis 


Routing  procedures  that  have  been  derived  to  optimize 
system  performance  in  the  sense  of  minimizing  the  total  delay 
accumulated  per  unit  time,  work  with  an  objective  function  of 
the  following  form: 


D, 


T 


L 

(i,k) 


^ik(^ik> 


(1.1) 


The  assumptions  inherent  in  (1.1)  are  discussed  by  Kleinrock 
in  [5].  denotes  the  average  delay/unit  time  on  link 

i-k  and  f^^^  denotes  the  total  link  flow  in  bits/sec.  An  extra 
term  may  be  added  to  include  the  effect  of  propagation  delays 
on  link  i-k.  These  previous  approaches  to  routing  have  employed 
closed  form  expressions  for  the  ® > derived  from  queue- 

ing theory  by  making  many  simplifying  assumptions.  This  approach 
to  the  routing  problem  has  been  taken  by  Kleinrock  [5],  Cantor 
and  Gerla  [2],  and  Schwartz  and  Cheung  [7j. 


The  departure  of  this  thesis  is  from  the  search  for  closed- 
form  expressions  to  finding  efficient  algorithms  to  estimate 
the  quantities  of  interest.  In  particular,  for  static  and  quasi- 
static routing  it  has  been  shown  in  the  previously  mentioned 
works,  as  well  as  others  like  Gallager  [3],  Agrew  [1],  and 
Segall  [6],  that  the  routing  procedure  should  be  based  on 
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knowledge  of  the  derivative  of  the  total  delay/unit  time  ^ 
of  messages  passing  through  a link  i-k  with  respect  to  the  total 


flow  Rather  than  differentiating  closed  from  expressions 

for  propose  processing  the  queues  at  the  links  to 

estimate  the  Dl,  (f.,  )'s  directly.  In  this  manner  we  can  dis- 
associate  the  optimality  of  a given  routing  procedure  from  all 
the  assumptions  necessary  to  the  closed  form  formulae  for  delay. 
In  this  thesis  w;  derive  three  different  estimation  procedures 
for  the  marginal  delays  ^)y  making  no  assumptions  as  to 

the  structure  of  the  queues.  However,  to  study  the  properties 
of  each  estimator  through  analysis  and  simulation  methods,  we 
make  very  specific  assumptions  about  the  structure  and  underlying 
statistics  of  queues  to  which  the  estimation  algorithms  are  to 
be  applied.  Hence,  in  the  following  paragraphs  we  motivate  the 
importance  of  directly  estimating  the  incremental  delays 
by  reviewing  the  previous  work  in  designing  routing  strategies, 
with  special  emphasis  on  those  results  relevant  to  quasi-static 
routing  procedures. 


1.3  Previous  Work 


The  most  common  model  for  routing  problems  in  data  networks 
is  that  derived  by  Kleinrock  [5].  He  makes  the  following 
assumptions : 


1)  FoisKon  arrivals  at  nodv.s 

2)  KxponenLial  distribution  of  mossaj^jo  lcMi^;tli 

3)  Independence  cf  arrival  processes  at 
different  nodes 

4)  The  "independence"  as.sum|ition  of  service 
times  at  successive  nodes.  Each  time  a 
message  arrives  at  a node  a new  service 
requirement  is  chosen  from  the  same 
exponential  distribution. 


On  the  basis  of  these  assumptions  he  derives  an  explicit  formula 
for  the  total  delay/unit  time  accumulated  on  the  (i-k)-th  link. 
If  denotes  the  amount  of  traffic  passing  over  the  (i-k)-th 

link  in  bits/sec.  and  is  the  capacity  of  link  (i-k)  in 

bits/sec.,  the  average  total  delay  will  be  given  by 


^ik  ■ ^ik 


(1.2) 


To  illustrate  how  Kleinrock's  result  in  (1.2)  is  used  in 
routing  problems,  we  outline  the  static  routing  scheme  presented 
by  Cantor  and  Gerla  [2].  The  problem  of  finding  an  optimal 
set  cf  routes  is  posed  as  a nonlinear  multi-commodity  flow  prob- 


lem, where  we  want  to  derive  the  flow  vector  f*,  whose  entries 


The  summation  in  (1.3)  is  taken  over  all  (i,k)  (lairs  that  are 
connected  and  y is  the  total  external  arrival  ri'ite  in  (jackets/ 
unit  time.  T is  interpreted  as  the  average  packet  delay.  The 
set  of  flows  that  satisfy  multi-commodity,  capacity,  and  non- 
negativity constraints  is  shown  to  be  a convex  polyredral  set 
and  hence  any  f in  that  set  may  be  expressed  as  a convex  combin- 
ation of  extremal  flows 

f = I L X.  1 (1.4) 

1 -1  3- 

1='  1=1 

Letting  vT(f*)  denote  the  gradient  of  the  objective  function  in 
(1.3)  evaluated  at  f*.  Cantor  and  Gerla  propose  an  algorithm 
that  finds  the  optimal  f = f-'’  in  the  sense  of  (1.3)  for  a given 

basis  of  external  flows  (cp^^^  ...  and  then  generates  a new 

(k+1)  (k+1) 

basis  vector  cp'  ' that  minimizes  (vT(f*)  .(P)  . (P^  is  the 

new  extremal  flow  that  will  help  us  reduce  T the  fastest.  The 

procedure  continues  iteratively  until  we  are  as  close  to  the 

optimal  flow  as  desired.  As  part  of  the  calculation  of  , 

we  generate  a set  of  routing  tables  that  realize  the  given  flow. 

Schwartz  and  Cheung  [7]  describe  a gradient  type  algorithm 

for  calculating  the  optimal  flow  vector  which  motivates  a possible 

stochastic  approximation  algorithm  for  quasi-static  routing. 

(m  n ) 

Let  f..  denote  the  total  bit  rate  on  link  (i-i)  and  f . . ’ denote 

ij  " ^ ij 

the  bit  rate  of  messages  with  source  m and  destination  n on  link 
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i-i.  Let  TI . denote  the  nropaioition  Lime  Lor  link  i-j  and  — 

ij  ■ p 

be  the  average  message  size  in  bits.  'Ihen  the  objective  function 
which  Schwartz  and  Cheung  define  is  the  average  message  dela,. 


T 


y 


T- 

(i,  j) 


f . . 

IJ 


(1.5) 


y is  the  expected  total  external  message  arrival s /unit  time  and 

, the  capacity  of  link  (i-j)  in  bits/sec.  If  NN  denotes  the 

number  of  nodes  in  the  network  and  v is  the  expected  number  of 

'mn 

arcivals/unit  time  at  node  m with  destination  n,  the  multi- 
commodity,  non-negativity,  and  capacity  constraints  on  the  flows 
are  stated  as  follows: 


NN 

L 

k=l 


(m,n) 

ik 


NN 

L 

f=l 


^(m,n) 
^ i 


(y  /m 

mn 

. 


i = m 

i^m,  i^n  (1.6) 
i = n 


f(m,n)  ^ Q 
ij 


(1.7) 


f.  . 

ij 


S 

(m,n) 


^(m,n) 

ij 


C.  . 
IJ 


(1.8) 


I 


Defining  comniodiLy  flow  veeLor  f,  whose  entries  are 

the  the  conservation  of  flow  constraints  (1.6)  may  be 

ij 

expressed  as 


A f = h . 

q 


(1.9) 


b is  a vector  whose  entries  are  either  0,  y /u,  or  -y  /u.  A 

mn  mn  q 

is  a matrix  consisting  of  submatrices  corresponding  to  each  (m,n) 


commodity. 


(1,2) 


(1,3) 


A = 

q 


(m,n) 


(1.10) 


* . ^(NN,NN-1) 


Given  a flow  satisfying  constraints  (1.6)  through  (1.8),  we 
can  obtain  a feasible  direction  of  descent  for  the  objective 
function  in  (1.5)  by  projecting  vT(f^)  onto  the  constraint  sur- 
face defined  in  (1.9).  Hence,  Schwartz  and  Cheung  propose  the 
iteration 


= f^  - hP  vT(c’-)  , 

q ^ » 


(1.11) 


where  h is  a step  size  and  a projection  operator  defined  by 


S' 


(1.12) 


r 


'P  T - 1 

P = 1 - A (A  A ) A 

q q q q q 


The  capacity  constraint  is  handled  implicitly  by  the  penalty 
function  in  the  objective  function  (1.5).  Schwartz  and  Cheung 
derive  an  h'  such  that  for  0 < h < h',  the  non-negativity  of 
flows  is  preserved.  The  actual  h is  determined  by  appealing  to 
the  convexity  of  the  objective  function. 

Now  we  recast  the  procedure  of  Schwartz  and  Cheung  [7]  to 
apply  to  a quasi-static  routing  situation.  Suppose  we  redefine 
the  objective  function  in  (1.5)  by  not  using  Kleinrock's  formula 


T = - L (D.  . (f.  .)  + f.  .mT:  .) 
y (i,j) 


(1.13) 


1. 


We  assume  that  the  current  vector  of  flows  is  f^  and  we  have 

5D 


available  estimates  for 


bf.  . 
ij 


and  the  Tl . are  known  constants 

£.fi 


Hence,  we  can  calculate  vT(f^)  and  apply  the  iteration 


-i+1 


= f^  - Tih'P  vT(f^)  , 


(1.14) 


where  h’  is  the  upper  bound  on  the  step  size  to  insure  non- 
negativity  of  flows  and  r;  is  some  scale  factor  0 ri  < 1 . 
liquation  (1.14)  could  bo  termed  a stochastic  approximation 
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aiguriLlun  since  our  values  ior  would  necessarily  be 

inexact  and  hence  our  gradient  vT(f^)  would  only  be  an  estimate. 

In  the  algorithm  offered  by  Cheung  and  Schwartz,  the  routing 
fractions  are  deti'rmined  by  knowledge  of  the  commodity  flows 
f^j'*n)  . In  addition,  the  procedure  is  centralized,  namely  we 
assume  the  routing  computations  are  performed  at  a special  node 
and  then  communicated  with  the  rest  of  the  network.  For  many 
reasons,  decentralized  algorithms  where  the  compu  ..ation  is  dis- 
tributed through  the  network  are  nu  e desirable.  We  next  discuss 
a quasi-static  routing  algorithm  de 'ived  by  Gallager  [3]  that 
not  only  is  decentralized,  but  works  directly  with  the  routing 
fractions . 

Gallager  uses  a static  model  with  stationary  inputs  and 

proposes  an  algorithm  that  seeks  to  minimize  the  total  delay 

in  the  network  specified  by  Eq . (1.1).  He  assumes  the  functions 

are  increasing  and  convex  U functions  of  the  flow  f • 

Let  t^(j)  denote  the  total  expected  traffic  at  node  i destined 

for  node  j and  ( j ) the  total  expected  external  arrivals  at 

node  i destined  for  node  j.  The  routing  variables  arc  defined 

as  (j),  the  fraction  of  traffic  t,(j)  that  is  routed  over 
i K i 

link  (i-k).  Conservation  of  flow  for  traffic  with  destination  j 
at  node  i is  expressed  with  these  variables  as 


(1.15) 


&D  bD 


I 

He  then  shows  that  a necessary  condition  for  i to  achieve  the  i 

i 

minimum  D_  is  \ 

T 1 


r 


Gal  lager's  algorithm  consists  of  two  parts;  a protocol 


between  nodes  to  calculate  marginal  delays 


dr^(j) 


rr  and  keep  track 


of  a number  of  sets  which  he  terms  and  a procedure  for 

up-dating  the  routing  variables  The  procedure  for  adjusting 
the  routing  variables  is  defined  as  a mapping  = A(«5)  that 
attempts  to  move  closer  to  the  optimal  equalibrii.m  condition 
specified  by  (1.19).  The  sets  B^(j)  denote  nodes  for  which 
^ik(j)  = 0 and  t’le  algorithm  is  not  permitted  to  Increase 
from  zero.  The  way  the  B^(j)'s  are  defined  insures  the 
"looplessness"  o.  routes  from  any  given  source  i to  destination  j, 
i.e.,  we  can  never  go  from  a node  i to  some  intermediate  node  m, 
back  to  node  i,  and  finally  to  our  destination  node  j. 

We  can  see  that  the  marginal  delays  fundamental 

to  Gallager's  procedure.  While  we  could  obtain  them  by  differen- 
tiating Kleinrock's  formula  t3],  it  would  then  be  necessary 
to  estimate  the  flow  f^j^.  Hence,  it  is  important  to  make  the 
algorithm  independent  of  Kelinrock's  assumptions  and  estimate 
directly  by  locally  processing  the  queue  for  traffic 
using  link  (i-k). 

We  next  review  Carson  Agnew's  discussion  in  [1]  of  the 
ARPA  scheme  for  routing,  since  his  analysis  reveals  the  reason 
for  the  sub-optimality  of  that  method  and  he  suggests  an  ARPA 
type  routing  strategy  employing  marginal  delays.  In  the  ARPA 
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procedure,  each  node  in  the  network  maintains  a table  whose  (i,j)-th 

entry  is  an  estimate  of  the  minimum  time  to  reach  the  j-th  node 

through  the  i-th  neighbor.  These  estimates  are  based  on  the 
queue  sizes  at  intermediate  nodes,  and  hence  the  time  it  takes  to 

empty  those  queues.  When  a message  arrives  addressed  to  the  j-th 

destination,  we  look  do^^m  column  j and  send  the  message  to  the 
neighbor  with  the  smallest  estimated  delay. 

Agnew  introduces  a simple  single-commodity  network  with 
input  flow  X to  be  split  into  n routes  with  X^  arrivals/sec. 
each  and  stationary  M/M/1  queues  (exponential  service  require- 
ments and  exponential  inter-arrival  times  all  mutually  independent). 
The  flow,  capacity,  and  non-negativity  constraints  together  with 
the  objective  function  corresponding  to  average  message  delay  are 
represented  as  follows: 

n 

X = L X.  (1.20) 

i=l  ^ 

0 < X^  < (1.21) 

T = f L X, (T,  + T!)  (1.22) 
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The  average  message  Length  is  s[>eeified  by  — . T\  denotes  the 
average  total  dclay/message  for  the  i ' th  queue  and  Tj  denotes 
some  remaining  constant  delay  to  get  to  the  destination,  such  as 
a propagation  time.  Agnew's  analysis  is  not  relevant  to  more 

general  networks  since  the  Tl  should  be  functions  of  X.. 

11 

For  this  simple  one  destination  model,  to  inpleraent  the 
ARPA  technique  we  would  have  a table  with  i ' th  entry  T|  + (l  + L^)gC^, 
where  denotes  the  number  of  messages  in  the  qu^ue  and  being 
serviced.  After  a long  time,  the  distribution  ot  traffic  would 
be  determined  by  the  equalibrium  conditions 


T.+t:=T.+t:  for  X.,A.>0 

1 1 J J 1 1 


T.  + t:  > T.  + T'.  for  X.  = 0,  X.  > 0 

3-  1-  J J 1 J 


(1.23) 


However,  these  do  not  correspond  to  the  conditions  for  obtaining 
a minimum  T in  (1.22),  which  Agnew  shows  to  be 


(X.T,)  + t:  = (X,T,)  + t:  X,,  X,  > 0 


bX.  1^  1 bX.  j i'  j 


X.  = 


0,  X . >0 

J 


(1.24) 
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Those  differing  equalibrinm  comiiLions  (1.23)  and  (1,24)  reflect 
the  difference  belween  system  and  user  u))  t imizat  ion . 

For  his  single-commodity  model,  Agnew  suggests  a way  to 
obtain  an  ARPA-like  scheme  that  will  approach  the  conditions  for 
system  optimality  defined  in  (1.24),  What  we  need  is  a quantity 


that  has  an  expectation  equal  to  the  marginal  delay  (X^T^). 

For  an  M/M/1  queue  he  shows  that  ,S.  = — ^ (1  3 b.)(l  + L./2) 

^ 1 1 1 

satisfies  the  desired  property.  Hence,  he  proposes  that  we  do 
ARPA-type  routing  with  revised  table  entries  + T|. 

In  [6]  Segal  1 proposes  for  a general  network  an  ARPA-like 
relating  strategy  that  uses  marginal  delays.  Suppose  that  the 
objective  function  we  wish  to  minimize  is  the  total  delay  accumu- 
lated per  unit  time  over  the  network  defined  in  (1.1).  If  we 
denote  by  r^(j)  the  average  bit  rate  of  external  arrivals  at  node 
i with  destination  j , and  given  a small  change  6r^(j)  in  r^(j), 
where  should  we  direct  the  extra  traffic?  Assuming  we  direct 
6r^(j)  on  a path  P from  i to  j , the  associated  change  in  the 
total  delay  is  given  up  to  first-order  terms  by 


60.^, 


(f,m)eP  -tm 


6r^(j) 


L D (f,  )6r.(j) 
,,  . -Cm"  -Cnh 

(t,m)cP 


(1.25) 
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Hence,  this  suggests  choosing  F so  as  1 1'  minimize 


routing 
scheme , 
list  the 
traffic 


procedure  motivated  by  (1.25)  is 
but  the  (i,j)-th  location  in  the 
estimate  of  the  minimum  ^ ^ 
to  destinaton  j through  neighbor 


6D 


1' 


Th  e 


analogous  to  the  ARPA- 
rouLing  table  would  now 
for  directing  extra 

i . 


1.4  Formulation  of  Thesis  as  Queueing  Theory  Pioblem 

In  all  of  1 lie  quasi-static  type  routing  algorithms  pre- 
sented, the  inc  'emental  delays  esse  itial  quantities. 

Rather  than  dif  crentiate  queuein^,  theoretic  forniulae,  with  all 
their  implied  siatistical  assumptions,  we  propose  estimating 

by  operating  on  the  record  of  the  queue  associated  with 
the  outgoing  link  from  node  i to  k.  We  are  interested  in  find- 
ing recursive  estimation  procedures  that  process  the  queueing 
record  and  converge  to  as  the  observation  interval 

becomes  sufficiently  large.  Hence,  our  problem  can  be  formulated 
in  the  context  of  single  server  queueing  theory.  In  our  case 
the  customers  are  identified  with  messages.  The  service  time 
becomes  the  transmission  time  for  the  message  due  to  the  finite 
capacity  of  the  communications  link. 

Hence,  in  this  thesis  we  derive  three  algorithms  that  pro- 
cess the  record  of  arrivals  and  departures  of  a queue  to  generate 
estimates  of  the  derivative  with  respect  to  arrival  rate  of  the 
average  total  delay  per  unit  time.  Since  we  make  no  assumption 
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as  Lo  Lhu  exact  runii  uL  the  (|ueiie,  we  ,ne  L li  1 1 ■ i i.: . L l'i.1  in  i ulai;  L 
techniques  that  are  as  insensitive  ti>  l .1 1 i ( i ca  ] assunpt  i ons 
of  the  queueing  process  as  possible.  llmavei',  the  per  lonaa'ic  1 ■ 
of  the  algorithms  can  be  analyzed  only  for  relatively  simple 
queues.  Consequently,  although  the  pro]:>i>,sed  algorithiiis  can  be 
applied  in  practice  in  very  g(?neral  situations,  their  explicit 
analysis  is  only  done  for  queues  like  M/M/1,  M/d/1,  etc. 

1 . 5 Summary  of  Thesis 

The  plan  of  the  thesis  is  to  examine  the  beliavior  of  tliree 
estimation  procedures,  which  we  term  the  customer-addition, 
customer-removal,  and  time-contraction  algorithms,  for  a variety 
of  queues.  The  algorithms  are  described  in  detail  in  Section  2.1 
but  we  may  say  now  that  each  procedure  corresponds  to  a different 
technique  for  imagining  a hypothetical  alteration  of  the  queueing 
record  to  reflect  a differential  change  in  arrival  rate  6X . In 
the  customer-addition  algorithm  we  conceptually  add  a customer 
at  a random  time  in  the  observation  period  to  simulate  an  increase 
in  arrival  rate.  In  the  customer-removal  algorithm  we  randomize 
the  conceptual  removal  of  customers  from  the  queue,  acliieving 
the  effect  of  a decrement  in  arrival  rate.  in  Llie  time-contraction 
procedure  we  redefine  the  arrival  times  of  customers  to  simulate  a 
compression  in  time  scale  and  hence  a differential  increasi'  in 
arrival  rate. 
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Sc’cLion  2 conLaiiis  the  main  I'e.sulLs  of  the  thesis.  First 


we  give  a detailed  description  of  the  way  tlie  notions  for 
altering  the  queueing  record  indicated  above  are  refined  into 
the  actual  estimation  procedures.  This  motivation  leads  then  to 
the  derivation  of  the  algorithms  and  their  reali/'ation  in  flow 
chart  form  is  ii  i cated.  The  analysis  of  the  algorithms  is 
performed  in  detail  for  special  queues  in  Sections  2.3,  2.6,  and 
2.8.  For  the  ci  ;tomer-addit ion  and  time-contraction  procedures 
we  are  able  to  p ove  asymptotic  unbiasedness  for  m M/d/1  queue. 
For  the  time-con  faction  and  customci. -removal  procedures  we 
define  the  calculation  of  the  asymptotic  bias  as  a power  series 
in  p,  the  utilization  factor  Xx,  in  the  case  of  M/g/1  queues. 
Employing  this  power  series  representation,  we  show  that  for  an 
M/D/1  queue  the  bias  for  the  customer-removal  algorithm  may  only 
contain  terms  of  third-order  or  higher  in  p.  We  also  show  that 
for  an  M/M/1  queue,  both  the  bias  for  the  time-contraction  and 
customer-removal  algorithm  contain  terms  with  powers  of  p of 
all  orders.  Since  the  calculation  of  the  variance  associated 
with  each  of  the  estimators  is  too  cumbersome,  we  derive  Cramer- 
Rao  bounds  for  each  algorithm  in  the  case  of  an  M/d/1  queue. 
Since  for  purposes  of  practical  implementation,  routing  calcula- 
tions are  secondary  to  the  actual  transmission  of  data,  it  is 
important  to  analyze  and  compare  the  storage  and  computation 
requirements  of  the  three  algorithms,  which  \-;e  accomplish  in  the 


final  ser.ion  of  Section  2. 
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In  SecLion  J we  present  the  results  ol  s imul  ,it  Lip',  .ill  this  . 


■i  - 


algorithms  for  an  M/l)/l  queue  and  the  cusLunier-renidv  il  and  tin., 
contraction  procedures  for  M/M/1,  D/M/I  , and  I'/M/l  queues.  '''  v- 
performance  measures  we  use  to  compare  the  algorithn^  are  the 
relative  bias  and  fractional  rms  error.  Sinci.  the  queueing 
record  is  segmented  into  busy  and  idle  (uriods  during  which  the 
server  is  occupied  and  unoccupied  respectively,  the  variable  we 
use  to  quantify  the  observation  interval  is  the  number  of  busy 
periods  N included  in  the  period.  Hence,  to  investigate  the 
convergence  of  the  algorithms,  for  each  queue  of  interest  we 
present  curves  of  the  fractional  rms  error  for  N = 10,  100,  and 
1,000  busy  periods  and  similarly  present  tables  of  the  fractional 
bias.  Employing  our  Cramer-Rao  bounds  for  the  M/D/1  case,  we 
find  that  all  three  algorithms  are  both  consistent  and  asymptot- 
ically efficient.  We  examine  the  robustness  of  the  customer- 
removal  and  time-contraction  algorithms  by  comparing  their  per- 
formance for  M/d/1,  M/M/1,  D/M/1,  and  U/M/1  queues.  The  only 
significant  difference  in  the  two  procedures  performance  occurs 
in  the  case  of  a D/M/1  queue,  where  the  customer-removal  pro- 
cedure does  worse. 
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SECTION  2 


THREE  ESTIMATION  ALGORITHMS 

2 . 1 Introduction 

The  main  goal  of  this  thesis  is  to  propose  and  evaluate 
algorithms  whith  process  the  record  of  a single  server  queueing 
system  to  estimate  the  derivative  of  the  total  delay/unit  time 
with  respect  tc  arrival  rate  X.  The  available  record  consists 
of  exact  knowlc  ge  of  arrivals  of  customers  to  i Te  queue  and 
their  departurto  after  service  is  completed.  Time  is  segmented 
into  alternate  intervals,  busy  periods,  during  which  the  server 
is  occupied,  and  idle  periods,  when  the  server  is  free.  The 
observation  interval  which  is  used  to  form  our  estimate  consists 
of  a number  of  busy  periods  and  the  intervening  idle  periods. 

A simple  thought-experiment  motivates  all  three  estimation 
algorithms.  Consider  our  single-server  queueing  system  with 
its  average  arrival  rate  of  X customers  per  unit  time.  For  a 
given  observation  period  T , if  we  can  compute  the  total  system 
time  S,  i.e.,  the  sum  of  all  customer's  service  and  waiting 
times,  then  the  average  delay/unit  time  is  given  by  D = S/T^^. 
Suppose  we  could  actually  alter  the  input  flow  by  some  6X . Then, 
on  the  basis  of  an  earlier  D,  by  computing  D-  for  the  next 
observation  period,  we  can  estimate  the  derivative  of  the  total 
delay/uni'  time  by  calculating 
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(2.1; 


However,  In  any  actual  queueing  system  it  would  be  undesirable 
to  cliange  flows  just  for  measurement  purposes.  Even  if  we 
could  implement  (2.1),  the  independent  statistical  fluctuations 
in  D and  D*  would  probably  make  it  a very  poor  estimator. 

Hence,  what  we  need  is  some  mathematical  formalism  for  an  imag- 
inary increment  in  flcv^;  6A,  which  will  allow  us  to  compute  the 
corresponding  change  in  delay  without  actually  perturbing  the 
arrival  rate. 

According  to  intuition,  an  increase  in  arrival  rate  should 

result  in  additional  customers  entering  the  system.  An  extra 

customer  arriving  in  a time  interval  T with  probability  c will 

h 

increase  the  effective  rate  by  6X  = extra  arrivals  are 

mutually  independent  events,  the  probability  of  two  or  more  cus- 
tomers will  be  of  second-order  in  e and  hence  of  second-order  in 
6X . Tlierefore,  only  the  effect  of  a single  extra  arrival  has  to 
be  considered  explicitly.  We  also  assume  that  the  arrival  time 
of  the  extra  customer  is  uniformly  distributed  over  the  observa- 
tion period  In  addition,  in  order  to  explicitly  compute  the 

change  in  total  system  time  over  the  observation  interval  due  to 
an  extra  arrival,  we  must  assume  that  the  additional  customer 
has  some  kno\'/n  service  requirement.  These  assumptions  allow  us 
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to  compute  iui  expected  increase  in  system  time  conditioned  on 
the  arrival  of  a new  customer,  and  the  resulting  estimation  pro- 
cedure will  be  called  the  customer -addition  algorithm. 

In  a second  algorithm,  an  incremental  decrease  in  the  effec- 
tive rate  X is  simulated.  lliis  is  done  by  assuming  that  each 
customer  arriving  to  the  system  is  allowed  to  indeed  enter  the 
queue  with  probability  1 -C,  and  is  eradicated  with  probability 
£,  independent!}  from  customer  to  customer.  In  this  way  we  simu- 
late an  arrival  process  with  rate  X(1  - e).  Hence  e is  determined 
as  follows : 


i - 


X(l-t)  = X + 


(2.2) 


c = 


X 


(2.3) 


We  then  estimate  X by  appealing  to  the  law  of  large  numbers.  If 


M'  is  the  total  number  of  customers  in  a period  T , then  XT  - M' 

h E 


Hence,  we  use 


-T, 


£ = 


M' 


6X 


(2.4) 


Again,  the  probability  of  removal  of  two  or  more  customers  from 


the  same  period  T is  second-order  in  6x  and  hence  the  reduction 

h 


of  total  system  time  that  has  to  be  considered  explicitly  is  due 
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to  removal  of  onl}’  one  customer. 


Tliir  reduction  cc'nsists  of  its 


own  system  time  and  the  effect  on  other  custc'iners.  llie  estima- 
tion procedure  motivated  by  this  second  techniiiue  for  making  a 
"virtual"  change  in  flow  bX  is  termed  the  cus tomer-removal 
algorithm. 

A second  characteristic  that  we  associate  with  an  increase 
in  arrival  rate,  besides  the  fact  that  more  customers  appear  in  a 
given  time  period,  is  that  there  is  less  time  between  successive 
arrivals  and  therefore  the  customers  are  more  "compressed" 
together.  To  make  this  argument  quantitative,  we  note  that  an 
average  arrival  rate  of  X customers  per  second  means  an  average 
inter-arrival  time  of  The  change  in  the  average  inter-arrival 

time  due  to  an  increment  in  X is  given  by 

6(t)  - - i (f)  . (2.5) 

If  T denotes  the  arrival  time  of  the  n-th  customer,  we  define  a 
n 

new  set  of  arrival  times  t'  = t (1  - 4^).  If  E(t  ,,  - r ) = ^, 

n n X n+i  n X 

then  - t')  = — (1  - • This  result  is  consistent  with 

n+1  n X X 

the  change  in  inter-arrival  time  predicted  by  (2.5)  due  to  an 
increase  in  flow  6X . We  compute  the  resulting  increment  in  system 
time  by  considering  first  the  fact  that  customers  arrive  a little 
earlier  and  second,  the  fact  that,  given  our  fixed  observation 

period  Tj.,  the  redefinition  of  arrival  times  results  in  tlie 
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trailing  edge  of  the  interval  being  contracted  and  leaving  ^ gap 
— T during  which  extra  customers  could  have  arrived.  The 
estimation  method  suggested  here  is  termed  the  t ime-contraction 
algorithm. 

The  present  section  contains  the  derivation  and  realization 
in  flow  chart  f.  rui  of  the  three  algorithms.  In  addition,  an 
extensive  analy^is  of  the  algorithms  is  performed.  We  give  a 
proof  of  the  as  nptotic  unbiasedness  of  the  custc  mer-addition  and 
time-contractio)  algorithms  for  a queue  with  Poisson  arrivals  and 
deterministic  service  requirements  (M/D/1).  The  asymptotic  bias 
behavior  of  the  customer- removal  and  time -con tract ion  algorithms 
for  M/G/1  systems  are  examined  as  a power  series  in  the  utiliza- 
tion factor  p = Xx.  For  an  M/D/1  queue  we  show  explicitly  that 
the  cus tomer-removal  algorithm  is  asymptotically  unbiased  up  to 
the  third  power  of  p and  give  a construction  to  prove  asymptotic 
unbiasedness  up  to  an  arbitrary  power.  Also  for  an  M/D/1  queue, 
Cramer-Rao  bounds  for  any  unbiased  estimator  of  the  delay  gradient 
are  derived.  They  will  be  used  in  Section  3 to  determine  the 
asymptotic  efficiency  of  the  algorithms  applied  to  M/D/1  queues. 

We  complete  the  present  section  by  analyzing  and  comparing  the 
storage  and  computation  reqi;irements  of  each  method. 
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2 , 2 l">ori  val  i I'n  and R^a  1 i / a M on  in  Flnv;  Pi  aj’ram  l-'orm  of 

Custoiner-Atl<.lIt:i()n  Al)’,oriLhm 


In  the  cus Loinor-acldLt: ion  algorithm  we  simulate  an  increase 
6X  in  the  arrival  rate.  The  following  assumptions  will  be  made: 

1)  The  probability  of  an  extra  arrival  in  the 

interval  T„  is  6XT,,. 

h t. 

2)  Each  extra  arrival  is  independent  of  all  other 
arrivals . 

3)  iTie  extra  arrival  is  uniformly  distributed 

over  the  interval  T^. 

E 

4)  The  service  requirement  of  the  extra  customer 
is  known;  we  denote  it  by  x. 

Let  Tj^  and  denote  the  duration  of  the  k-th  busy  and  idle 
periods,  respectively.  Let  6S  denote  the  increase  in  system  time 
over  N busy  periods  associated  with  the  arrival  of  an  extra 
customer.  We  let  6S  denote  the  expected  increase  in  system  time 
associated  with  an  increase  in  arrival  rate  6X  and  conditioned  on 
the  record  of  arrivals  and  departures.  By  conditioning  on  the 
random  arrival  time  t being  in  each  Tj^  and  we  can  compute  6S 
as 

6S  = E(6S iQueueing  Record)T^6X 


6S 


< L E(6S'|tcT.  , Queueing  Record) 
(k=l 


E(6S'|telj^,  Queueing  Record) 


N-1 

Z 

k=l 


(2.6) 
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The  T 6A  outside  the  brockets  is  tin-  probability  of  an  extra 

arrival.  The  increment  in  system  time  due  to  an  increase  in  flow 

5A  is  zero  il  no  additional  arrival  occurs.  The  factors  Tj^/T^^ 

and  I,  /T  represent  the  probabilities  of  t being  in  the  k~th  busy 
K h* 

and  idle  periods,  respectively.  This  is  a consequence  of  the 

assumption  that  t is  uniformly  distributed  over  T . Since  we  are 

h 

interested  in  tl^e  derivative  of  the  total  delay/unit  time  with 
respect  to  the  flow  rate,  our  estimator  is  given  by 


6' 


(2.7) 


We  focus  next  on  the  calculation  of  E(6slte:Tj^,  Queueing 
Record)  and  E(6SjtcIj^,  Queueing  Record).  These  expected  incre- 
ments in  system  time  are  composed  of  the  average  effect  on 
existing  customers  plus  the  average  system  time  of  the  additional 
customer.  In  considering  additional  arrivals  in  a busy  period, 
we  can  distinguish  between  effects  on  the  customers  in  that  busy 
period  and  interactions  with  succeeding  busy  periods.  First  we 
examine  the  part  of  E(6S  \ tcTj^,  Queueing  Record) , call  it  ASj^, 
that  comes  from  considering  the  k-th  busy  period  in  isolation. 

To  facilitate  discussion,  the  following  notation  is  defined. 


T.  = Arrival  Liiuu  oi’  i-t.Ii  cusLumcr  in  tliv  busy  period 
(relative  to  the  start  of  the  busy  interval) 

X.  .^Service  reciuirement  of  the  i-lh  customer 
1 — ' 

S.  .^System  time  of  the  i-th  customer 

S*  ^ System  time  the  i-th  customer  would  havt  had  if 
an  additional  customer  arrived  at  time  c 

M ^Number  of  customers  in  the  busy  period 

A ^ System  time  of  the  additional  customer 

t ^Arrival  time  of  the  additional  customer 

T ^Duration  of  the  busy  period. 


(2.8) 


We  break  up  the  calculation  of  the  expected  increase  in  system 
time  into  expectations  conditioned  on  an  arrival  in  the  interval 
[t.,  t.,,]  for  i = 1 . . . M.  T,  is  zero  and  t.,,,  is  defined  as 
the  duration  of  the  busy  period  T.  Then 


•4 


AS  = E (L  S*  + A - L S lte[0,T]  ] 
n n ' 


M 

= L E [S*  + A - L S^lte[T.,T._^^]l  Pr  [tf  ltc[0,T]] 

i=l 

(2.9) 
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By  the  assumption  of  uniformly  dis tril)utecl  arrival  time 


Pr  [te[T^,T^_^^]  lte[0,T]  ] 


T . , - T . 

1+1  1 


(2.10) 


The  system  time  of  a given  customer  is  equal  to  his  service  plus 
waiting  time.  Tlu  waiting  time  is  equal  to  the  sum  of  the  ser- 
vice requirements  of  those  who  entered  the  busy  period  before  him 
minus  his  arriva  . time.  Hence,  the  system  time  of  the  n-th  cus- 
tomer is  given  by 


n 

S = S X.  - T 
n 1 n 

1=1 


(2.11) 


Now  consider  the  new  total  system  time  due  to  an  arrival  at  time 

tUTi.T.^P' 


M 

L S*  + A 
n=l 


1 

S S + 
n 


M 

S 


n=l 


(S  +x)  + 
n=i+l  \ n=l 


1 

S X 


n 


+ X - t . 


(2.12) 


The  first  term  represents  the  first  i customers  whose  system  times 

are  unaffected  by  the  new  arrival.  The  second  term  shows  that 

each  customer  ahead  of  the  new  arrival  will  suffer  an  additional 

delay  x.  The  final  term  represents  the  system  time  A of  the  new 

customer.  Since  conditioned  on  being  in  the  random 

variable  t is  uniformly  distributed  over  that  interval,  taking 

appropriate  conditional  expectations  in  (2.12)  yields 
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K [i:  s,;  t A - 


S [ T . , T 
n ' 1 


i+1 


i 

X + L 
n=l 


X 

n 


+ (M-i)x  . (2.13) 


Substituting  (2.13)  and  (2.10)  into  (2.9)  results  in  the  follow- 
ing expression  for  AS: 


AS 


= X + 


M 

i=l 


L 

n=l 


X 

n 


M 

+ !:  (M-i)x 

i=l 


(2.14)  ] 


The  first  and  second  terms  represent  the  service  requirement  of 
the  new  customer  and  his  expected  waiting  time,  respectively. 

The  third  term  is  the  expected  delay  suffered  by  the  existing 
customers  due  to  the  new  arrival. 

For  the  special  case  where  all  the  service  requirements  are 
the  same,  namely  x^  = x,  Eq.  (2.14)  simplifies.  This  can  be  seen 
by  direct  substitution  of  x^  = x into  Eq . (2.14),  but  it  will  be 
illuminating  to  re-examine  the  computation  that  led  to  (2.14). 
Suppose  our  extra  customer  arrives  at  time  t after  the  k- th  and 
before  the  (k+l)-st  customer.  The  M-k  customers  ahead  of  the 
new  arrival  suffer  an  additional  delay  x.  According  to  the  rule 
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described  by  (2.11),  tlie  sysLem  time  of  the  new  customer  is  ^jiven 
by  (k+l)x-t.  Hence,  for  an  arrival  at  time  t,  there  is  an  addi- 
tional system  time  Mx-t+x.  This  amount  of  Lime  is  equivalent  to 
allowing  the  extra  customer  to  wait  till  the  end  of  the  busy 
period  and  then  i>e  served.  Since  t is  uniformly  distributed  over 
an  interval  Mx,  v;e  have  t = Mx/2.  Therefore,  AS  is  given  by 


AS  =7,-  Mx+x  , 


(2.15) 


We  complete  the  calculation  oi  E(6S|tcTj^,  Queueing  Record) 
by  looking  at  the  additional  system  time  that  may  result  from 
one  busy  period  overlapping  onto  another,  ITie  following  will 
hold  for  arbitrary  service  times  x^.  No  matter  where  an  addi- 
tional customer  arrives  in  the  k-th  busy  period,  that  period  will 


be  extended' by  the  extra  service  time  x.  The  value  of  x relative 
to  the  following  idle  period  durations  will  determine  the  number 
of  succeeding  busy  periods  that  will  be  affected  by  an  arrival  in 
Tj^.  However,  the  number  is  always  finite.  If  x < no  follow- 


ing busy  periods  suffer  additional  delay.  If  < x £ ’ 
only  the  (k+l)-st  busy  period  is  affected.  Tlie  exact  effect  on  a 
given  busy  period  j in  the  future,  depends  on  how  much  an  arrival 


in  T,  causes  the  (j-l)~st  busy  period  to  overlap  onto  the  j~th  busy 
K 

period.  F'or  example,  if  x > then  each  customer  in  will 


suffer  an  additional  delay  (x-Ij^)  . 


Letting  denote  the  number 
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of  customers  served  in  the  k- th  busy  period,  the  preceding  reason- 
ing leads  to  the  following  rule  for  computing  E(6S  |teTj^,  Queueing 
Record) ; 


E(6S 

Queueing  Record) 


AS, 


'''  \+l^^  “ ^k  - ^ - ^k  ^k+1 


- 


f+1 


j-1 


L •141 


AS,  4-  L X - S I,  ^ K + I-  ip+A 

j.l  ^+j\  „.o  '"■^/j-0  J-0 


(2.16) 


The  last  relation  in  (2.16)  refers  to  the  case  when  an  arrival  in 
Tj^  affects  (-64-1)  busy  periods  into  the  future. 


To  complete  our  description  of  the  customer-addition  algor- 
ithm, we  must  now  evaluate  the  average  increase  in  system  time 
E(6S  jtelj^,  Queueing  Record)  associated  with  arrivals  in  idle 
periods.  The  effect  of  an  arrival  in  on  the  j -th  busy  period 
again  depends  on  how  much  the  (j-l)-stbusy  period  slides  onto  the 
j-th  busy  interval.  As  the  effect  on  customers  in  T^  is  not 
independent  of  the  exact  arrival  time,  wc  must  average  over  all 
times  t in  the  k-th  idle  period.  Let  denote  the  time  instant 


corresponding  to  the  beginning  of  the  (kfl)-st  busy  period.  If 


X < 1,  an  arrival  in  1,  aflecLs  only  the  (k+l)-sL  busy  period. 
— k“T'i  k. 

Hence  E(6Slt€l,  , Queueing  Record)  is  computed  as 
K. 


Queueing  Record) 


= X + 


max 


/"k+l  *'*  ■*' '■  "“k+1^  1, 


for  X < I,  , 1 . 
— k+1 


(2.17) 


Here,  x + t - o,  ,,  is  the  amount  that  the  additional  customers 
service  time  overlaps  onto  the  (k+l)-st  busy  period.  The  "max" 
is  necessary  in  the  lower  limit  of  integration  since  our  arrival  t 
must  be  in  the  interval  Ij^.  If  x < then  - x represents 

the  earliest  time  at  which  an  arrival  ran  occur  and  influence 
the  (k+l)-st  busy  period.  By  a simple  change  of  variable, 
t'  = t - the  dependence  on  disappears.  The  relation 

(2.17)  may  be  generalized  to 


E(6Sltelj^,  .0  ^ 

n • =x  + L / M,  (x  + t'  - S. ) ^ 

Queueing  ■ ^ f-I  S -xl  ^ h 

Record)  ^ ^ ^ 


for  S I,  , . <’h<  < L I,  , . . 

j-i 


(2.18) 


Equation  (2,18)  corresponds  to  the  case  when  an  arrival  in  1 
influences  -l+l  succeeding  busy  periods.  The  value  of  the  inte- 
grals in  the  summation  are  given  by 


(x  - X I,  -S.)  -1,  S . -X 
2 k j'  k j 

(x-S  . ) ^ S . -X  > -1, 

J J k 

(2.20) 

Employing  (2.18)  and  (2.16)  in  (2.6)  and  (2.7)  we  can  con- 
' ceive  of  a processor  which  up-dates  an  estimate  for  the  de  1 ay  gradient 

at  the  end  of  each  busy  period.  Let  denote  the  time  from  the 
start  of  the  observation  period  to  the  end  of  the  k-th  busy  period. 
Let  denote  the  estimate  for  the  delay  gradient  based  on 

k+1  busy  periods.  Let  be  the  incremental  expected  delay 

suffered  by  the  (k+l)-st  busy  period  due  to  an  additional  arrival 
in  the  current  queueing  record.  Hence,  the  u[i-dating  at  the  end 
of  the(k+l)-st  busy  period  assumes  the  form 
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S.)  ^ 


(2.21) 


D) 


k 


(k+1 ) 


^(k)  ^ "^k  ' 


Since  x is  finite,  we  look  back  in  time  a finite  number  of  busy 
and  idle  periods  to  compute  Suppose  that  the  busy  and  idle 

periods  are  numbered  consecutively  from  the  beginning  of  the 
observation  time  and  let  i^  denote  the  index  of  the  most  recent 
idle  period.  Let  n^  denote  the  index  of  the  first  idle  period  at 
which  an  arrive!  with  service  requirement  x can  influence  the 
current  (ii  i 1)  st,  busy  period,  L >nce,  we  need  only  store  idle 
period  and  busy  period  information  for  (I^  ...  _l_^)  and 

(I  . . . T.  At  the  end  of  every  idle  period,  n^  must  be 

n^+l  i-r+i  1 


up-dated  to  reflect  how  far  back  we  must  look  to  compute  effects 
on  the  newest  busy  period.  Additional  simplification  is  possible 
due  to  the  fact  that  both  E(6SlteIj^,  Queueing  Record)  and 
E(6SlteI,  , Queueing  Record)  include  an  x term.  Since  the  incre- 
mental  expected  system  time  due  to  an  arrival  in  either  a busy  or 
idle  period  is  weighted  in  the  estimator  by  the  probability  for 
arrival  in  that  time  slot,  we  can  add  x at  the  end.  These  ideas 
are  realized  in  the  flow  diagrams  of  Figures  2.2  and  2.3. 

Figure  2.1  pictures  the  relationship  between  n^  and  i^,  and 
defines  the  variables  that  appear  in  the  flow  diagrams  of  Figures 
2.2  and  2.3. 
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r 


L 


Hi  + l 


4 h 

I. 

1, 


T. 


Time 


^1+1  I 
I 

_ J 


ij+l  = Number  of  the  current  busy  period. 

Hj.  = Number  of  the  first  idle  period  at  which  an 
arrival  with  service  requirement  x can  cause 
customers  in  busy  period  to  suffer 

additional  delay. 

l,l'  = Variables  denoting  the  elapsed  time  from  the 
beginning  of  the  observation  period  to  the 
end  of  busy  periods  i^,  respectively. 


A,  S,  C = Auxilliary  variables. 

q = Index  of  current  busy  period 

= Arrival  time  of  n-th  customer  in  i-th  busy 
period  relative  to  the  beginning  of  that 
busy  period. 

x^  = Service  requirement  of  n-th  customer  in  the 
i-th  busy  period. 


= Number  of  customers  in  the  i-th  busy  period. 

/s 

D'  = Delay  gradient  estimator. 

/s 

= Delay  gradient  estimator  minus  service 
requirement  of  additional  customer  x. 


Figure  2.1  Queueing  Record  Structure  and  Definition  of  Variables 

Relevant  to  Flow  Chart  Realization  of  Customer-Addition 
Algorithm 
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Wc  now  pose  the  question  of  whether  the  eustoiner  addition 
algorithm  can  be  generalized  to  be  applicable  to  a wider  class 
of  queues  than  those  where  all  customers  have  the  same  service 
requirement  x.  As  formulated,  the  algorithm  is  limited  by  an 
assumption  of  a fixed  service  requirement  x for  the  additional 
customer.  Hence,  we  can  conceive  of  extending  the  algorithm  by 
doing  a final  averaging  over  x, 


D' 


j 


00 

/ B(x)  If  (x) 
x=0 


dx 


(2.22) 


^ s 

B(x)  denotes  the  service  time  density  and  (x)  refers  to  the 
unnorraalized  incremental  delay  as  a function  of  the  assumed 
extra  customer  service  requirement  x.  While  possible  in  principle 
the  scheme  implied  by  (2.22)  is  unacceptable  for  practical  reasons 
The  evaluation  of  (2.22)  necessitates  saving  the  entire  queueing 
record  and  doing  all  our  processing  at  the  end. 


Only  when  the  service  time  density  consists  of  a discrete  set 
of  values  would  it  be  reasonable  to  implement  (2.22).  In  this 
situation  B(x)  is  given  as  a train  of  impulses. 


L 

B(x)  = L p,  6(x-x,  ) (2.23) 

k=l 


Then  (2.22)  would  become 


r 


■wnw 


D'  = 


E k=l 


6S  . V 

f“k  6T  <’‘k> 


(2.24) 


Hence,  for  each  Xj^,  k = 1 ...  L v\?e  would  process  the  queueing 
record  in  parallel,  employing  the  algorithm  given  in  the  flow 
diagram  of  Figure  2.2. 

2 . 3 Proof  of  Asymptotic  Unbiasedness  of  Customer-Addition 
Algorithm  for  an  M/P/1  Queue 

We  now  examine  the  bias  of  the  customer-addition  algorithm 
as  the  number  of  busy  periods  in  the  observation  period,  N, 
become  unbounded.  For  the  algorithm  to  be  asymptotically 
unbiased  we  must  prove  that 


lim 

N-*«> 


Jl-  M U = ^ 


6X 


bX 


(2.25) 


i. 

« - 


where  D is  the  average  total  delay/unit  time. 


bD 


To  check  (2.25)  we  must  first  define  the  quantity  ^ . The 
average  total  delay/unit  time  D is  equal  to  X times  the  average 
total  delay/customer  . Hence,  ^ may  be  expressed  in  terms  of 


5D 
c 

bX 


as 


V ^ ■ 


T,  b I) 

^ = D + X — 
bX  c bX 


(2.26) 
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Thu  average  total  clclay/custoiner  is  expressible  in  terms  of  the 
average  service  time  x and  the  average  waiting  time  w as 


D=x+w.  (2.27) 

c 


Hence 

time 


5D 


can  be 


’ DX 
and  servic  e 


reformulated 
recjuirement . 


in 


terms  of  the  average  waiting 


= '■  + w + A ^ (2.28) 

C OX 

We  can  evaluate  the  above  expression  for  all  queues  for  which  an 
explicit  form  of  the  waiting  time  distribution  is  available. 

Reviewing  the  assumptions  inherent  in  the  customer-addition 
algorithm,  we  can  expect  that  the  procedure  will  be  asymptotically 
unbiased  in  the  case  of  an  M/D/1  queue.  The  descriptor  "M/d/1" 
means  the  arrival  process  is  Poisson,  and  the  service  require- 
ments deterministic.  Since  all  customers  in  an  M/D/1  queue  have 
the  same  service  requirement,  the  assumption  that  the  additional 
customer  has  a fixed  service  time  is  harmless.  The  two  other 
assumptions,  uniform  arrival  time  distribution  for  the  extra 
customer  and  the  probability  density  for  the  arrival  of  extra 
customers,  are  both  consistent  with  a Poisson  arrival  process. 

For  a Poisson  process  the  probability  density  for  the  time  of 

occurrence  of  the  i-th  event  given  k > i events  did  occur  in  [0,T] 
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I 


is  uniform  on  the  interval.  Let  p(k,X)  denote  the  probability 
of  k arrivals  in  an  interval  T given  that  the  arrival  rate  is  X. 

If  we  let  (6XT)  be  the  probability  that  k additional  customers 
arrive  in  an  interval  T due  to  an  increase  in  rate  6X , p(k,X+6X) 
must  satisfy 

p(k,X+6X)  = p(k-i,X)(6XT)^  + p(k,X)|l  - (2.29) 

After  some  manipulation,  dividing  both  sides  by  6X  and  taking 
the  limit  as  6X  approaches  zero,  we  obtain 


bp(k.l) 

bX 


Tip(k-l,X) 


p(k,X)]  . 


(2.30) 


By  direct  substitution  we  can  verify  that  the  Poisson  process 
formula  for  the  probability  of  occurrence  of  k events  in  a time  T 
given  below  satisfies  (2.30), 


p(k,X) 


(XT)^i~^^ 

k! 


(2.31) 


Motivated  by  the  preceding  arguments,  we  proceed  to  prove 
that  the  customer-addition  algorithm  is  asymptotically  unbiased 
for  an  M/D/1  queue.  Since  for  an  M/d/1  queue,  the  average  wait- 
ing time  is  given  in  [4]  as 
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i 


■ V 


(2.32) 


2U-P)  ’ 


where  p = Xx  is  the  utilization  factor,  formula  (2.28)  dictates 
that  we  must  prove  the  expectation  of  our  estimate  (2.7)  converges 
as  N " to 


(2.33) 


NM 

Since  by  the  lav  of  large  numbers  we  have  lim  T = lim  — , 

N-.od  ^ N-.oo 

interchanging  tie  limit  and  expectation  operations  in  (2.25)  we 
must  show  that 


...  / 1 „ 6S  \ _ OD 

N->co  nM  / 


(2.34) 


Employing  (2.6),  we  can  break  the  problem  of  computing  E 
into  evaluating  the  expectation  of  two  types  of  terms  as  below. 


E 77  = L E[E(6S\teT,  , Queueing  Record)!,  ] 


+ S E(E(6S|tel,  , Queueing  Record)!,  ] 
k=l  ^ 


(2.35) 


The  terms  E(6S|teTj^,  Queueing  Record)  and  E(6S|t€lj^,  Queueing 
Record)  are  given  by  Eqs . (2,16)  and  (2.18),  respectively.  Since 
the  queueing  record  is  given  by  the  number  of  customers  served  in 
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each  busy  period  (Mj^  . . . and  tlie  idle  period  durations 

(I^  ...  ^),  we  further  break  the  calculations  by  first  condi- 
tioning on  (M^  . . . and  averaging  over  (1^^  ...  and  tV.en 

averaging  over  . . . Mj^)  as  below. 

E[E(6S  jteTj^,  Queueing  Record)Tj^] 

= E[E[E(6S\t€Tj^.Mj^  . . . ‘ ^ (2-36) 


E{E(6S \telj^,  Queueing  Record)!^^] 

= ElE[E(6S\tcIj^,M^  --Vl  •••  ^N-l^^k^^l  ' ' ' ^ (2-37) 


We  organize  our  calculations  by  first  computing  6S  and  6S 

1 ^1 

defined  below  as 


and 


A E[E(6SltcT^,M^  . . . 


(2.38) 


6Sj  ^ E[E(6S  itfIj^,M^  . 


M 1 
N’  1 


and  L'lU  ii  ^v  ucra]  idling  our  ruaulL  Lo  5S,,  and  6S  specified 

k k 

similarly.  Finally,  we  sum  6S.  and  6S,  over  all  k's  and 

k k 

average  over  (M^  . . . . 

Calculation  of  6Sj 


We  calculate  6C  by  first  illustrating  the  thinking  involved 
when  the  number  of  busy  periods  included  in  the  observation  per- 
iod, N,  is  thre  • and  then  generalizing  the  procedure.  Figure  2.4 
depicts  the  que  eing  record  for  N=3  and  a partiL  oning  of  (1^,12) 
space  into  thre  regions:  R2,  and  R^.  According  to  (2.16), 

T^E(6SlteT^,  Queueing  Record)  is  given  by 


r 


TiE(6Sltcq,  -(4S,I,  + (x-ipM^T, 


X < or 


(q.ipcRi 


< X < I2  or 


(I^,l2)eR2 


X > + I2  or 


^ + (x-I  j^)M2Tj^ 

+ (x-I^-l2)M3T^  (I^,l2)eR3 


(2.40) 


The  key  to  computing  the  desired  expectation  is  in  noting  which 
region  of  the  (1^,12)  space  corresponds  to  each  term. 
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avei'ngcd  over  the  whole  space.  is  averaged  over 

£ X and  ( x-I  ^ - 1 )M2T  is  averaged  over  +1^  £ x.  Sub- 
stituting = M^x  and  employing  (2.15)  for  AS,,  we  break  up  the 
expectation  by  averaging  each  term  over  tlie  appropriate  region. 

6S,^,  = -f  M^x^')  + Mj^M^xECx-Ij^  jx  > I^)  Pr  (x  > ) 

+ 1^  lx  > d l2>  • ^2^  (2.41) 


We  can  now  ext  ul  the  arguments  that  led  to  (2.  <f)  to  the  case  of 
an  arbitrary  n uber  of  busy  perie  s N.  The  necessary  condition 
that  there  be  ii  contribution  to  tl  e incremental  delay  due  to  the 
effect  of  an  arri\'al  in  on  the  (k+l)-st  bus)  period  is 


k 

L 


I . 
J 


X . 


(2.42) 


Hence,  '2.41)  generalizes  in  the  case  of  N busy  periods  to 

2 122  ^ 

6S  - (M  X + - M^x  ) + L M R xE  X - LI. 

1 k=l  d ^ j 

( I 

Pr  L I.  £ X . (2.43) 

(j=l  I 

We  now  proceed  to  defire  the  statistics  of  quantities  which 
we  need  to  calculate  (2.43).  For  any  M/g/1  queue,  tlie  idle 
period  lengths  are  independent,  identically  distributed  expon- 
ential random  variables  with  parameter  A [4].  Hence,  their  sum 


k 

L 

j=l 


I. 

J 


< X 
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I'.as  a gaiTima  distribution.  Tlie  density  and  distribution  function 
associated  with  the  sum  Y of  k idle  period  durations  are  defined 


below. 


■i  « 


Y = L I. 

J = 1 ' 


(2. 44) 


^v(y) 


k-1  k^-Xy 


y X e 


(k-1)! 


(2.45) 


Fay)  - 1 - 

^ j=0 


(2.46) 


To  calculate  (2.43),  we  need  to  evaluate  E(Y\y  < x) . 


E(Yly  < x)  = 


0 


yfy(y)  ^y  - j 1 


-Xx  ^ jX^ 


Pr  [y  < x] 


i I 


1 - 


-Xx  (Xx)^ 

e Z.  — 


(2.47) 


j=0 


J! 


Using  (2.47)  and  (2.46),  Eq.  (2.43)  becomes 


N-1 


6S^  “ (M^x^  + I mJx^)  + L I 


k=l 


-Xx  (Xx) ^ 

c L 


j=0 


I 


+ e 


, ,k-l  k 
-Xx  X X 


(k-D!  J 


(2.48) 
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Tlie  first:  term  repretiOuLs  the  efleet  of  the  additional  arrival  on 
the  first  busy  period  and  the  terms  in  tlic  summation  show  effects 
on  the  remaining  N-1  busy  periods. 

Calculation  ofc  6S 
1 

The  preced  ,,  g procedure,  of  oxaraining  each  type  of  term 

separately  and  .rnposing  conditions  on  the  space  or  (I^  ...  1^  j) 

such  that  the  t rm  appears  in  E(6S  |tf  i'p  Queueing  Record)!^,  may 

be  applied  to  c aputing  6S  . A more  compact  st  rement  of 

^1 

Eqs.  (2.18)  thr  ugh  (2.20)  for  E(6h  Queueing  Record)  will 

make  the  identification  of  the  relevant  terms  clearer. 


r t+i 

X + L 


E(6S 


Queueingv 

Record 


+ (l-vp 


!ki 


(x-Sj ) 


2( 


-0 

L 


j=l 


k+j 


< X 


t+l 

L 

j=l 


k+j 


V 


for 


(2.49) 


From  (2.44),  there  are  three  terms  in  E(6SiteI^,  Queueing  Record)! 
to  consider. 


xl. 


(2.52) 


^1+j 


S.)I, 


(2.53) 


M 


(2.54) 


The  first  term  (2.52)  always  appears  and  hence  is  averaged  over 


the  whole  space  of  (1^^  ...  We  next  consider  the  j = l terms 


specified  by  Eqs.  (2.53)  and  (2.54).  Iliere  is  no  condition  on 
(1^  ...  Ijyj  needed  to  guarantee  contributions  to  tlie  incremental 
delay  due  to  the  effect  of  an  arrival  in  1,  on  the  second  busy 


period.  Relation  (2.51)  imiilios  that  expression  (2.53)  will 
appear  i i'  Ij^  < x and  expression  t2.54)  if  j > 1. 

.i 

L I < X is  the  necessary  condition  that  arrivals  in  1,  affect 
1=2 

the  (jH)~st  busy  period.  Taken  together  with  (2.501  and  (2.51) 
this  implies  th*  following  rule  for  computing  tha  expectation  of 
terms  (2.53)  an’  (2.54)  for  j > 1. 

1 j ( j ) 

Average  M.,  . (x  - y I,  - L I,  )li  Over  Z j , < x (2 . 55) 

^ 1=1  - ( 1=1  ^ ‘ 

j 

L 1.  > x 
1=1 


j 

L I,  < X 
1=2 

This  discussion  of  the  conditions  for  the  appearance  of  all 
the  terms  in  E(6Sjtelj^,  Queueing  Record)!^  is  summarized  in  the 
formulation  of  its  expectation. 


(2.56) 


M 


Average 


i+i 


(X  - 


J 

L 

1=2 


^1>' 


Over 


xEI. 


N-1 

Z 


Mi^.E[(x-2 


S.)It  is.  + I,  <x1  Pr  [S.  ■+!.  <xl 
1 ' J 1 ' J 1 


N-1 
+ L 
j=l 


^ E[  (x-S  .l''lS.+I.>x,  S.  <x]Pr[S.+lT>x,  S.<x 
2 J ' J 1 J J 1 J 


(2.57) 
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*r  • • --K- 


V ♦ 

' 


s.  = 

J 


J 
L 
f =2 


j = l 


.58) 


is  fhe  unconiK 1 1 onal  mean  of  '.n  exponential  random  variable 
with  parameter  > and  'uence  i'  computation  for  the  j-1 

A 

terms  in  tlu  se.nu  ijti  ms  are  luniped  U'gctlu  r and  the  result  listed 
b<.  > <jv. . 


(2.59) 


i > 


■jlv's  term  i.s  tha  'part  of  the  ! m r <i  li.en  ta  I sys';ej;i  time  oue  to  the 
effect  on  tlie  busy  period  following  The  terms  for  j > 1 

represent  con  t ri '.'utions  d'.tr.  to  effect.s  c".  busy  periods  more  tlian 
one  res'.oved  froiti  To  com()ule  tlie  term?  in  (2.57)  for  j > 1 

we  can  rewrite  the  exp(;ctations  implied  in  (2.55)  and  (2.56)  as 


Average  M 


■14  ) Pb 


j \ 1 ,1  I j 

1-1  ’')b  "2 


/ — 


I 


(2.60) 
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•'=1 


Averaj;;t*  1,  + ) Uvc-r 


(2.61) 


I 


By  defining  two  random  variables  we  can  reformulate 

the  evaluation  of  the  expectations  implied  by  (2.60)  and  (2.61). 


Y,  = I 


1 n 


(2.62) 


Y.  = 


- 


I.  = 


J 

£ I, 


t=l 


(2.63) 


The  joint  density  for  Y^^,  is  computed  from  the  density  for  the 
su.m  of  j-1  independent  exponential  variates. 


-xy.  /'xJ-hy,-y.)J-V<>'2-y 


(j-2)! 


^Y^,Y2^^1’--'2^  (j-2)!  ^ 


-2  “^^’2 


y2  2:  Y|  1: 


(2.64) 


We  nov.'  define  two  regions  in  (Yj,Y2)  space, 


) -1  ‘ 


UY,  ,Y.):  Y, 


(2.65) 


r 


IKj) 


Yi  1 


0 y, 


7?tKT 


1 , , -AX.  —X X — X X ^ X 

- (1  - f ) - xe  - Xc  -^:pjyy 


(2.70) 


. XX  Y 

(Y^Ri)  = I ; Yi  dy^dy^ 

0 y^ 


2 — Xx. 

-o  (1-e  ) 

X 


2 — Xx  2— Xx 

~ xe  - X e 
A 

Pr  ;Rj) 


,-xkJ;2  2x'-x'--^3 

""  Jo 


(2.71) 


’ r ,.  ,T  ^ 1 > •*-  ^ 


(Y  y,m,)  - ,'  f V,v„  ‘l’'2  ‘ , , 

1 2 ' 1 ^ J 12  — ,R^1  — tiy2<iyi 


( 2 — Xx.  2 — Xx  2— Xx 

7 -X  (1  - e ) - - xe  - x e 

I ^ 


/I  Ax^  X X 1 

+ — (-  (l-O  ) - xc  1 


-■y-  (,)H1)  Iv  -^TTQTr 


___ 


61 


(2.72) 


Pr  = r 


>'2=-'^  yi=y2“^ 


(j-1)! 


(Yi1R2>  = I 


V ^Y,Y^^1’^2^ 

^2  1’  2 


y2=-->--  yi=y2-^ 


^1  Pr  [R2]  ^y\^y2 


(X  + p 


Pr  tR2] 


(Y^  IR2) 


^Yi, Y2^^1’^2^ 


Y2=^  Yj^=y2"^ 


^1 MR2T 


liiiii:!  ,,2^1 


0^)T"  ■*  X \ (j^+j)le^’‘ 

X 

Pr  1R2] 


(2.75) 


' * »: 


(Y21R2) 


Y2-^  yi"y2~^ 


2 

Pr  R 1 
2 


dyjdy2 


(Ax)^  ^ / I l.-Ax 

= ilZDUlLlKil.- 


PrlR^'t 


(2.76) 


(Y2IR2)  = 


Y2=x  y^^y2  ^ 


fy  V ^Yl 
2 1’  2 

'2  PriR^] 


dy;^dy2 


(Ax)-^"^  ,2^2  2 ,--AX 

F^iR^l 


(2.77) 


(Y1Y2IR2) 


yz=x  y^=y2-x 


YiY2 


V (YijYo^ 


pFOg — “>'i‘''’2 


(,2  + ii±Ll  , + J X)^xx 

pTTrP 


(2.78) 


Substituting  (2.69)  - (2.78)  into  (2.67)  and  (2.68)  and  using 
(2.59),  we  can  finally  evali ate  the  expeclation  of 
E(5SUcI^,  Queueing  Record)  Ij  over  (Lj  ...  l.^_j)  outlined  in 
(2.57). 
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+ ± 4 Jjjlil  + ill  2-Xx 

i 2'  2 ^ X -•-  2 e 

A A 


N-l 

^ i,  ”-j 


j-2 

xe  L 


1+2 


j-2 


1 _^— Xx  -'^‘'  ( X X ) ' " ^ J_  — Xx  ^ (Xx) 


X-  ,Io  (^+2)!  ^2-  (t+3)! 


M I 1^2L)—  21  I 


(2.79) 


Calculation  of  6S„  , 6S  and  E 


To  use  our  results  (2.48)  and  (2.79)  for  6S„  and  6S 

n ^1 

respectively,  in  order  to  derive  the  mean  of  the  unnormalized 

estimator  formulated  in  (2.34),  we  note  that  in  computing  6S_ 

^k 

and  6S_  , only  the  N-k  idle  periods  following  T,  enter  into  the 
^k 

averaging.  Hence,  we  use  our  answers  for  6S^  and  6S*  adjusted 

1 ^1 

to  correspond  to  a N-k+1  busy  period  case.  By  the  preceding 
argument,  the  expectation  of  the  unnormalized  estimator  over 
dl  • • • is  given  by 


'^1 
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I (Mj^x  4 j t L 

k=l  j=l 


L 

k=l 


M.M 

J 


j+k\ 


<f + hv+i<'<  + ^ 


N-2 

I 

k=l 


N-k 


L 

j=2 


b. 

J 


N-2  N-k 


The  behavior  of  and  are  made  clearer  tiy  replacing  each 

^ X 

summation  in  their  definition  by  e minus  some  quantity. 


1 


= X (x 


k V X X 


T=k 


(Ax)- 

-t! 


. 2— Ax 

+ X e 


k-1 


t'k-1) ! 


(Ax) 


t+2 


u = ^ 21  T 

j “ ^ I A (t+2)! 


(>-x) 


T+3 


.2  t.j-i  <«)!  ( 


(2.84) 


Evaluation  of  Limit  in  (2.34) 


Having  almost  evaluated  the  expectation  of  the  unnormalized 
estimator,  we  are  nearly  ready  to  examine  the  limit  in  (2.3^). 


We 


6S 


complete  the  expectation  of  by  averaging  over  (M^  . . . M^^)  . 

— 2 ~ 

This  amounts  to  replacing  by  M,  IL  by  M , and  noting  that  due 

—2 

to  the  independence  of  the  M's,  = M . Carrying  out  the 

expectation  over  the  M's,  changing  the  order  of  summation  in  the 
double  sums,  and  dividing  by  NM  we  are  left  with 
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*:  *7:'  . r.> 


1 6S  2 , 1 2 - L (N-k)a 

— E ^ = X + 2 ^ k=l  ^ 

NM  M 


^ NziljL^l  (,  - l))[ 

^ IxM  ^ ^ ) 


+ i L (N 
N j=2 


N-1 

■j)b  + ^ L (N-j)C  . 

J ^ j = 2 -' 


examine  the  behavior  ot  (2.85)  as  N - « we  mus'-  evaluate 


the  following  limits: 


1-  1 N-1 

lim  1 ^ , 

— T,  ka, 
N-”  N , . k 

k=l 


T , N-1 

i Z jb. 

N-‘->  N ^"2  J 


1 N-1 

i I ](.. 

N"'“  N j_2  j 


(2.89) 


->  *:-ir 


(2.90; 


N-1 
L C. 
j=2 


(2.91) 


To  prove  that  the  three  limits  (2.86)  - (2.88)  are  zero,  it  is 
sufficient  to  show  that  the  unnornialized  infinite  sums  are  finite. 
The  infinite  sums  implied  by  (2.86)  - (2.91)  are  evaluated  by 
switching  the  order  of  summation  and  looking  for  terms  that  cor- 
respond to  the  exponential  power  series.  The  results  are  expressed 
in  terms  of  p = Xx. 


_ , 2.1  , 1 2, 

L ka.  = x [-  p + -7  p 1 

k=l  ^ 


(2.92) 


S jb,  = 


x^e"^  2 3 0 

= 2 — ^3  b e 


- h p^  + p^  + 2e^  - 2pe^  - 2} 


(2.93) 


I jC,  = 


2 -p 
X e 


[pc^  - + 1 - 


(2.94) 


\ = 2 


(2.95) 


CO 


b. 

J 


(2.96) 


j=2 


2 ^ 


2 -p 
X f 


X 

X 


(c 


-P 


1) 


(2.97) 


Examination  of  the  power  series  in  p that  corresponds  to  (2.92)  - 
(2.94)  shows  that  each  is  a bounded  function  of  p on  [0,1],  We 
are  interested  in  p on  [0,1]  since  the  statistics  of  the  queueing 
system  are  stationary  for  this  range.  The  boundedness  of  (2.92)  - 
(2.94)  implies  that  the  limits  (2.86)  - (2.88)  must  be  zero. 

To  complete  the  description  of  (2.83),  we  list  from  [4] 
expressions  for  the  first  and  second  moment  of  the  number  served 
in  a busy  period  for  an  M/D/1  system. 


” ■ 1-p 

(2.98) 

g2.  >- 

/ -X  _ \ ^ ^ P 

(2.99) 

Employing  our  knowledge  of  the  limits  (2.86)  - (2.91)  provided  by 
(2.92)  - (2.97)  and  using  (2.98)  and  (2.99),  we  can  calculate  the 
limit  suggested  by  (2.85). 
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2 


(2.  lUO) 


I ini 
N 


i.  r ^ ^ 

- '■■  6X  ( ? • (l-p) 


+ f (l-p)  + ^ X^J^) 


Multiplying  Eq . (2,85)  Ey  X and  with  some  minor  rearranging  we 
obtain 


2(l-p) 


(2.101) 


2(l-p)' 


This  is  the  desired  delay  gradient  for  an  M/D/1  queue  derived  in 
(2.33).  Hence,  on  the  basis  of  the  thinking  leading  to  (2.34), 
we  have  proven  the  asymptotic  unbiasedness  of  the  customer- 
addition  algorithm. 

2.4 

Since  the  calculation  of  the  exact  variance  associated  with 
the  customer-addition  algorithm  is  too  cumbersome,  we  derive  a 
Cramer-Rao  bound.  If  we  have  an  observation  vector  R,  a para- 
meter A we  want  to  estimate,  and  a conditional  density  (R|A), 

the  Cramer-Rao  bound  for  the  variance  of  any  unbiased  estimator 
S(R)  of  A is  stated  as  follows: 


Cramer-Rao  Bound  for  Customer-Addition  Algorithm  in 
Case  of  M/d/1  Queue  ' 


In  the  case  of  an  M/D/1  queue,  the  observation  vector  which 
the  customer-addition  algorithm  works  with  is  the  concatenation  of 
two  sets  of  variables.  If  Y denotes  the  total  observation  vector. 


then  Y ^ (Yj^'.Y^)  where  Y^^  consists  of  (M^  ...  M^)  and  of 
(I^  ...  ^).  We  know  that  is  independent  of  and  1^^ 

independent  of  1 for  all  i j,  k ^ -t.  The  length  of  idle  period 
is  determined  by  an  "end"  effect  in  the  k-th  busy  period. 

Hence,  the  only  conceivable  place  we  could  find  statistical  depen- 
dence is  between  and  Ij^.  We  resolve  this  question  by  consider- 
ing the  density  for  the  length  of  idle  period  conditioned  on 
= m^.  The  dynamics  of  successive  waiting  times  in  a queue  are 
described  by  the  recursion 

cu  . 1 = max  {0,  u:  + x ~ 6 j with  initial  condition  oo,  =0. 

n+1  ’ n n n ^ 

(2.103) 

X is  the  n-th  service  requirement  and  6 the  inter-arrival  time 
between  the  n-th  and  (n+1 ) -s t customer . If  there  are  ra^  customers 
in  the  first  busy  period,  the  following  relations  must  hold. 
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+ Xj  - 0.  > 0 


i = 1 ...  rtij^  - 1 


(2.104) 


uj  + X - 6 <0 

nij^  nij 


(2.105) 


I.  = 0 - (x  + u;  ) 

i rrij^  rtij^  nij^ 


(2.106) 


For  any  M/g/1  queue,  the  unconditional  density  for  6^  is  exponen- 
tial with  parameter  A.  Hence,  the  density  for  conditioned  on 
customers  in  the  first  busy  period  is  related  by  (2.105)  and 


m 


(2.106)  to  the  density  for  0^  , conditioned  on  0^  being  greater 
than  the  sum  of  the  m^^-th  service  and  waiting  time. 


-i  T 


P(t) 

6 1 6 > X 4 w 

m^  m^  — m^  m^ 


-X  (x  -f-  w ) 
m^^  m^^ 


Xe 


-X  (t  - (x  -*-0:  ) ) 

1 "^1 


for  T > X -f-  u: 

— m^^  m^ 


(2.107) 


From  (2.106)  and  (2.107)  we  calculate  the  conditional  density  of 
as 


-XI. 


P(ll) 

hlMi 


- e 


(2.108) 


Hence,  since  the  density  for  conditioned  on  mj^  customers  being 

served  in  the  first  busy  period  is  identical  to  the  unconditional 
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density,  Mj  and  1 are  s tat  is  t it  a 1 1 y independent  random  variables. 
This  indepcmlence  is  a proi)orty  of  the  'memoryless  ' inter-arrival 
density.  Based  on  the  preceding,  argunieii  t s , all  the  variables  in 
the  observation  vector  Y ^ f'M^^  . . . 1 j ...  mutually 

independent . 

Hence,  the  joint  density  of  Y may  be  expressed  as 
N N-1  -XI. 

P(Y)  = IT  Pr  [M.  - m.}  TT  Ae  ^ . (2.109) 

i=l  ' j = l 

For  M/G/1  systems,  queueing  theory  has  calculated  the  probability 
of  k customers  being  served  in  a busy  period  [4] . The  result 
for  an  M/D/1  queue  is  given  in  [4]  as 


Pr  (M.  = m.  ] 
1 1 


ra.  - I 

(m^p)  -m^p 


m . 


(2.110) 


For  the  moraent,  we  pretend  that  the  parameter  of  interest  is  p 


and  using  (2.110)  we  reunrite  (2.109)  in  its  terras. 

N / N 


N ra . 

P(Y’p)  = I TT 

i=l 


m^  - 1 


L ra , - N - 

i=l  ^ 
p e 


N-1 

- ^ L I 


L m.  p 

1-1  ^1  =*  i-i 

X 


.1 


(2.111) 


Wf  U'L  y denote  tlie  delay  grailieiit  For  p mi  [0,1]  Fq . (2.3J) 

specifies  a 1-1  correspondence  between  y and  p.  Relation  (2.33) 
may  be  inverted  to  find  p as  a function  of  y. 


p = 1 -VI  - 


(2y  - x) 


We  can  evaluate  the  second  partial  derivative  of  the  logar- 
ithm of  the  joint  density  required  in  (2.97)  by  applying  the 
chain-rule  of  differentiation 


tn 


2 

6 t-n 


Performing  the  above  manipulations  and  employing  the  two  follow- 
ing expectations: 


E L m.  = 

i-i  ^ i-" 


N-l 

E E I.  = 


We  can  evaluate  (2.92)  for  the  customer-addition  algorithm. 


Var  (y  - y)  > 


(2.116) 


2 2 

2<_J2 

(N  - 1 +p)  (1  - p)^ 


The  result  (2.116)  behaves  as  expected  for  p near  zero  and 
one.  As  the  utilization  factor  p nears  1,  the  queue  becomes  non- 
stationary. The  means  and  variances  of  variables  such  as  the 
number  in  the  system,  the  waiting  time,  and  the  number  served  in 
a busy  period  become  infinite.  Hence,  any  estimation  algorithm 
which  is  a function  of  these  queueing  variables  might  be  expected 
to  diverge  as  p goes  to  1 . If  we  conceive  of  o approaching  zero 
by  fixing  x and  letting  X go  to  zero,  the  average  idle  period 
duration  becomes  unbounded.  In  addition,  Var  goes  to  zero 

as  p -•  0.  Hence,  since  the  queueing  variables  become  "known"  as 
p -»  0,  it  is  reasonable  to  expect  that  the  varian  of  the 
estimator  goes  to  zero  for  p = 0. 


2 . 5 Derivation  and  Realization  in  Flow  Diagram  Form  of 
Customer-Removal  Algorithm 


Since  the  function  U(x)  for  the  average  delay/unit  time 
accumulated  by  the  queue  is  continuously  differentiable  on 

1 n . , 

0<X*^  — (0  < p < l)  , the  imi  t defining  ^ 3 given  X'"^ 

X 

is  independent  of  the  direction  from  which  X approaches  X*. 


"X 


O(^)  - t)(X^v) 
X - X* 


x=x* 


liiii 

X - X* 
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(2.117) 


In  the  customer-addition  algoriLiim  we  let  X = X*  + 6x  and  allov/ 

6X  to  approach  i>ero  tlirough  positive  values.  For  the  customer- 
removal  algorithm  we  equivalently  let  6X  go  to  zero  through  nega- 
tive values.  We  simulate  a decrement  in  arrival  r\ce  6X  by 
removing  customers  from  the  queue  with  probability  L and  computing 
the  resulting  decrement  in  total  system  time. 

The  value  for  c is  motivated  by  the  fact  that  the  expected 

change  6XTg  in  the  number  of  customers  arriving  in  an  interval 

T„  caused  by  a decrement  6X  in  incoming  flow  equals  the  negative 
h 

of  the  expected  number  of  customers  removed  by  Bernoulli  trials, 

M*C.  Here  M'  is  the  total  number  of  customers  arriving  in  a per- 

(1) 

iod  1'  . Hence,  c = - — ; — . Letting  6S . denote  the  change  in 
h M j 

system  time  of  the  i-th  busy  period  due  to  the  removal  of  the 
j -th  customer,  the  expected  change  in  system  time  6S  in  a time 
Tg  due  to  a decrement  in  flow  6X  and  conditioned  on  the  queueing 
record  is  formulated  as 


( 


D' 


1 E(6S iQueuein^  Record) 

T 6X 

E 


_1_ 

M' 


N 

V 

i-1 


M. 

1 


j = l 


(2.119) 


We  compute  65^^^  by  working  with  more  microscopic  quantities. 

Let  c'^  . denote  the  amount  of  system  time  saved  for  the  M-th 
m,  1 

customer  in  the  i-th  busy  period  by  the  removal  of  the  n-th  cus- 
tomer in  that  busy  period.  Since  the  removal  of  the  n-th  cus- 
tomer can  have  no  effect  on  customers  that  preceded  him, 

. = 0 for  m = 1 ...  n-1.  Hence,  can  be  computed  as 

m,i  J 


6s: 

J 


(i)  _ 


M. 

1 

L 

m=j 


m,  1 


(2.120)  I 

i 


We  now  develop  a systematic  procedure  for  calculating  the 
c”  .'s.  To  simplify  the  notation,  we  drop  the  i denoting  the 

K , 1 

index  of  the  busv  period.  Let  oj  , S , and  x denote  the  waiting 

time,  system  time  and  service  requirement,  respectively  of  the 

n-th  customer  in  the  busy  period.  Let  d and  a denote  the  cor- 

responding  departure  and  arrival  time  of  the  n-th  customer. 

Since  the  system  time  the  n-th  customer  saves  by  the  removal  of 

the  n-th  custcmier  is  S , we  have  = S . In  considering  the 

n n n 

effect  of  rc'moving  the  n-th  customer  on  the  (n  + l)-st  customer, 
either  a new  busy  period  begins  with  the  (n  + l)-st  customer  or 
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V A 


•>  - 


the  (n  + l)-st  customer  remains  part  of  the  busy  period  formed  by 

customers  1 to  n-1.  The  conditii'n  for  customer  n+1  beginning  a 

new  busy  period  is  that  the  arrival  time  of  the  (n  + 1) -sc 

customer  is  greater  than  the  departure  time  d^  ^ of  the  (n-l)-st 

customer.  In  this  case,  customer  n + 1 will  save  its  waiting  time 

GO,,.  Ifd  customer  n + 1 does  not  start  a new  busy 

n+1  n-1  n+1’  ^ 

period,  and  saves  an  amount  of  time  since  it  need  no  longer 

wait  for  customer  n to  be  served.  This  rule  for  , , is  summar- 

n+1 

ized  by  the  following; 


for 


for 


a , , > d 
n+1  — n-1 


d 1 > a , , 
n-1  n+1 


(2.121) 


Relation  (2.121)  is  more  succinctly  stated  as 


,n 


■'n+l 


= d 


max 


u 


n+1’ 


Vi’ 


(2.122) 


Noting  that  max  [a  , , , d ,]  = -min  (-a  , -d  , , may  be 

^ n+1  n-1-^  n+1’  n-1  ’ n+1  ^ 

restated  in  a final  form  as 


"n+1 


= min  I GO 


n+1 


(2.123) 


Similar  reasoning  to  that  employed  in  calculating  applies 

to  the  computation  of  c/\  The  roiaoval  of  customer  n either  causes 

customers  m and  m-1  to  be  in  the  same  busy  period,  or  customer  m 

may  begin  a new  busy  period.  The  removal  of  customer  n causes 

customer  m-1  to  save  system  tim.e  Hence,  customer  m-1  departs 

at  an  earlier  time  d , - . . If  this  new  departure  time  for 

m-1  m-1 

customer  m-1  is  greater  than  the  arrival  time  a of  customer  m, 
customers  m and  m-1  remain  in  the  same  busy  interval  and  customer 
m is  saved  a system  time  . . However,  if  a > d , - , , 

customer  m begins  a new  busy  period  and  saves  its  waiting  time 
These  relationships  are  summarized  in  the  following  rule  for  com- 
puting c" : 


..n 

■'m 


c”  , for  d , - > a 

m-1  m-1  m-1  m 


d 1 - a = o)  for  a > d ^ - c”  , 
m-1  mm  m m-1  m-1 


(2.124) 


This  rule  may  be  expressed  more  compactly  as 


C - min  [C  - , a:  j 
m m-1’  m’ 


(2.125) 


Hence,  the  algorithm  for  computing  the  may  be  summarized 


as 
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c"  = s 

n n 


c’^ , , = min  tx  , w . } 
n+1  t»  n+1 


(2.126) 


c'^  = min  [C^  T , to  ^ m - n + 2 ...  M 
m m-1’  ra' 


M is  the  number  of  customers  served  in  the  given  busy  period. 

The  customer-removal  algorithm  is  completely  specified  by  (2.126) 
and  the  following  form  for  the  delay-gradient  estimator  derived  by 
substituting  (2.120)  into  (2.119). 


D'  = 


M’ 


N M 
Z L 
i=l  n=l 


M 

E 

ms=n 


_n 


1 


N 


m,  1 


M * ^ 

^ i=i 


M 

L 

m=l 


m 

E 

n=l 


,n 


(2.127) 


The  second  form  suggests  calculating  and  summing  for  n = 1 . . . m 
when  the  m-th  customer  arrives.  Hence,  to  calculate  the  inner  two 
summations  in  (2.127),  we  need  only  M variables  to  store  as 
j is  varied.  This  idea  is  realized  in  a flow  diagram  for  the 
customer-removal  algorithm  in  Fig.  2.6.  The  variables  in  the  flow 
. : 4tt  are  <iefined  in  Fig.  2.5. 
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D > 


' = Current  estimate  for  delay  gradient 

M'  = Current  total  number  of  customers  in  observation 
period 


TS  = Running  sum  of  service  times  in  most  recent  busy 
period 


X = Service  requirement  of  most  recent  customer 

T “ Arrival  time  of  most  recent  customer  relative  to 
beginning  of  busy  period 

w = Waiting  time  of  most  recent  customer 


j = Index  of  most  recent  customer  in  current  busy 
period 


M = Total  number  of  customers  in  most  recent  busy 
period 


= Storage  location  for  C.  as  j is  varied 


M M 

S = L L C for  niost  recent  busy  period.  It  is 
n=l  m=n 

the  cumulative  system  time  saved  by  the  removal 
of  each  customer  in  the;  current  busy  period. 


Figure  2.5  Definition  of  Variables  for  Customer-Removal  Algorithm 
Flow  Diagram 
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and  Initial-  4.-*-  x 

ize  4.  LJ ^ 


Compute  and  ___  i_-[ 

Add  Ccntri-  . 
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i=l  . 


J S—  S + 


j-1 

^ M+M' 
M-»-  M’+M 


M'+M 


■“  Update  Delay 

Gradient  Estimate 
S at  End  of 
Busy  Period 


Figure  2.6  Flow  Diagram  for  Custoiior  Removal  Algorithm 


2 . 6 Calculation  of  Asymptotic  Bias  foi  Customer-Removal 

Algorithm  for  M/G/1  Queues 

We  now  investigate  the  asymptotic  properties  of  the  customer- 
removal  algorithm  by  first  interpreting  the  terms  in  the  estimator. 
For  a given  busy  period,  the  inner  two  summations  in  (2.127)  may 
be  grouped  into  two  terms  representing  the  sum  of  all  the  service 
times  of  the  customers  in  that  busy  period  and  the  cumulative 
service  time  saved  by  all  other  customers  due  to  the  removal  of 
each  customer  separately.  If  denotes  the  system  time  of  the 

j -th  customer  in  the  i-th  busy  period,  the  customer-removal  delay 
gradient  estimator  may  be  expressed  as  follows: 


D' 


N 

Z 

i=l 


M. 

1 

L S 

j=l 


(i) 


N 

z 


+ 


N 

S 

i=l 


P. 

1 


M. 


N 

Z 

i=l 


M. 

1 


where  F.  is  defined  by 

1 ■' 


(2.128) 


We  examine  the  asymptotic  beliavior  of  the  mean  of  the  esti- 
mator specified  in  (2.128)  by  interchanging  the  expectation  and 
limit  operation.  By  appeal!. tg  to  the  law  of  large  numbers,  the 
limiting  form  of  the  estimator  as  N becomes  unbounded  is 

Urn  ED'  = E lim  6'  = D +1:  (2.130) 

N ® N ® M 

where  is  the  average  system  time  per  customer  and  P is  the 
expectation  of  the  quantity  defined  in  (2.124).  For  an  M/G/1 
queue,  the  mean  of  is  independent  of  i since  the  ' s depend 
on  waiting  times  and  service  times  which  are  statistically 
independent  from  one  busy  period  to  anotlier.  M denotes  the 
average  number  of  customers  served  per  busy  period. 


Using  the  general  relation  derived  in  (2.26)  that  the  delay 
gradient  is  equal  to  plus  X we  can  formulate  the  asymptotic 

bias  of  the  customer-removal  algorithm  as 


b = lim  ED ' 

H ..  cr 


bX  M 


X 


bD 
c 

bX 


(2.131) 


We  can  break  up  the  calculation  of  P by  conditioning  on  M = i for 
i = 2 . . . ® and  find 
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where 


P = L Q(i)l\ 
i = 2 


(2.132) 


f ^ Probability  that  i customers  are  served  in 

i * 

a busy  period 


^ n 

Q(i)  t Z T,  E(C^lM=i) 
n=l  m=n+l 


(2.123) 


(2.134) 


For  an  M/G/1  queue,  D is  given  in  [4]  as 


^ I ^ 2(1  -P  ) ” ’ 


(2.135) 


where 


(2.136) 


denotes  the  variance  and  x is  the  mean  of  the  service  time 
b bD 

distribution.  Hence,  X -7'-  may  be  expressed  as 


^ TT  - ' 


(1  + 


2 2(1  -p)^ 


(1-p)  f (1  + O 


(2.137) 


> .V.JI 


The  average  number  oi  customer!-  served  per  busy  period  is  given 
in  [4]  by 


M = 


1 -p 


C2.138) 


The  z-transform  for  the  probability  density  of  the  number  of 
customers  served  per  busy  period  is  described  in  [2.4]  by  the 
following  functional  equation: 


F(z)  = zB*[X  - XF(z)j  . 


(2.139) 


1^ 


B*  is  the  one-sided  Laplace  transform  of  the  service  time  density 
and  F(z)  is  defined  by 


F(z)  = L f z". 
n=l 


(2.140) 


Having  described  the  quantities  that  compose  the  asymptotic  bias 
b,  one  may  express  it  as 


b = (1  -p) 


L Q(n)f  - 
n=2  " 


^ .CH12)  (k±lj  ^k^l 

k=0  ^ 


(2.141) 


Now  suppose  P in  (2.132)  may  be  expressed  as  a power  series  in  p. 
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•■V 


P = L Q(n)f 
n=2 


* - 

j=0  J 


(2.142) 


then  based  on  (2.141),  to  prove  the  customer-removal  algorithm  is 
asymptotically  unbiased,  we  must  show  that 


for  n = 0,1,2,  ...  (2.143) 


In  addition,  using  (2.142),  the  power  series  for  b is  given  by 


= [ttQ  - 2 x(l+C^^)lp  + L (a  -a 

j = l 


I x(l  +C^)(j+l)]pj'^\ 

(2.144) 


For  an  M/D/1  queue,  we  were  able  to  prove  (2.143)  for 
n = C,  1,  2,  which  shows  t^iat  for  this  case,  the  estimation  algor- 
ithm is  asymptotically  unbiased  at  least  up  to  third  order  in  p. 
Although  we  could  not  prove  (2.143)  for  arbitrary  n,  based  on  our 
intuition  we  believe  that  for  M/d/1  the  algorithm  is  asymptotically 
unbiased.  We  have  also  examined  (2.143)  for  M/M/1  queues.  For 
this  case  it  turns  out  that  (2.143)  does  not  hold  even  for  n = 0, 
which  shows  that  the  asymptotic  bias  contains  terms  of  order  p. 

In  the  remainder  of  this  section  we  give  the  detailed  proofs 
indicated  above  for  M/D/1  and  M/M/1  queues.  The  main  part  is 
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the  calculation  of  Q(i)  defineci  in  (2.134),  so  that  we  divide  the 


proofs  into  scvetal  steps: 

a)  Since  by  (2.121),  is  given  as  a function  of  waiting 

n+k 

times  and  service  rejquireinents 


C = min 
n+k 


lx 


n 


to 


n+1  ’ ^^n+2 


oc 


h+k 


(2.145) 


we  calculate  the  joint  density  of  (x^,  •••  <J^p+k^  conditioned 

on  M. 

b)  Evaluate 

c)  Proof  of  (2.143)  for  M/d/1  for  n = 0,  1,  2. 

d)  Calculation  of  first  and  second  order  terms  in  p in  the 
asymptotic  bias  for  M/M/1  queues. 

Calculation  of  Joint  Density  of  •••  ‘^n+k^ 

Conditioned  on  M 

We  approach  the  problem  of  deriving  p(x^,o;^_l_^  ... 
by  deriving  the  joint  density  of  the  service  and  inter-arrival 
times  conditioned  on  M customers  being  served  in  the  busy  period. 
Let  denote  tlie  inter-arrival  time  between  the  j-th  and  (j+l)-st 
customer.  By  the  rules  for  conditional  probabilities,  we  can 
break  up  the  joint  density  of  (Xj^,  ...»  Xj^,  ...  0^^)  con- 

ditioned on  M as  follows: 


p(Xj^  . . . Xj^,  0J^  . . . = 


p(MlXj  ...  x^,e^  ...  0j^)p(xj  . . . Xj^,0^  ...  0j^) 

Pr  [M] 


(2.146) 


p(M\x^  ...  •••  0j.j)  is  either  one  or  zero  depending  on  whether 

the  variables  (x^^  ...Xj^,0j^  ...  0j^^)  satisfy  the  constraints  such 
that  exactly  M customers  are  served  in  the  busy  period.  With 
no  conditioning,  the  service  times  and  inter-arrival  times  are 
mutually  independent  random  variables.  Since  our  queue  is 
M/G/1,  the  inter-arrival  times  are  exponential  random  variables 
with  parameter  X.  Hence,  p(x^  ...  Xj^,0^  ...  is  given  by 


M M -xe. 

p(x,  ...  x^.e,  .••.6m)  " ^ B(x.)  ir  Xe  (2.147) 

1 1 ” i=l  ^ i-1 

The  conditions  on  the  x ^ ' s and  0 ^ ' s that  guarantee  M customers 
are  served  in  the  busy  period  come  from  requiring  that  the  waiting 
times  of  customers  2 ...  M are  greater  than  zero  and  satisfying 
a terminal  condition  that  customer  M+1  falls  outside  the  busy 
period.  Karlier  wc  defined  the  waiting  time  of  the  k-th 
customer  as  the  sum  of  the  service  requirements  of  preceding  cus- 
tomers minus  his  arrival  time  relative  to  the  start  of  the  busy 
period.  This  relative  arrival  time  may  be  expressed  as  the  sum 
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of  the  first  k-1  in  L <>r-tir  r Lv.i  1 Limes  0 . . Hence,  Lhe  condition  for 

J 

the  waiting  timi!S  of  cnsttMiiers  2 ...  M being  greater  than  zero  may 
be  expressed  as 


k-1 


X . 

J 


k-1 

V 

j=l 


e . > 0 

J - 


for  k 


2 ...  M . 


(2.148) 


For  the  (M+l)-st  customer  to  fall  outside  the  busy  period,  the 
arrival  time  of  the  (M+l)-st  customer  relative  to  the  start  of  the 
busy  period  must  be  greater  than  the  departure  time  of  the  M-th 
customer.  Noting  that  the  M-th  customer  departs  when  all  M ser- 
vice requirements  have  been  satisfied,  we  define  a dummy  variable 
to  state  the  terminal  condition. 


u^^l  = I X - Z e < 0 (2.149) 

j=l  j=l 


Hence,  the  joint  density  of  M service  times  and  inter-arrival 
times  conditioned  on  M customers  being  served  in  the  busy  period 
is  as  follows: 


K=1 


®1  " ' ■ ®M 


M M -X9. 

■n  B(x.)  IT  Xe 
i=l  ^ i=l 


k 

M-1  0 < e,.  i s X. 

^ *1  J 

j=i 


k-1 


'M 


M 

> L 

j = l 


M-1 
- Z 

j = l 


0.  . 

J 


X . 
.1 


(2.150) 


As  defined  by  (2. 140),  denot  cs  the  probability  that  a busy 
period  has  M customers. 

We  can  now  calculate  p(x^  ...  •••  W working  with 

the  linear  relationship  between  (6^^  ...  0j^)  and  ... 

implied  by  Eqs.  (2.148)  and  (2.149).  The  inverse  relations  for 
the  6j  ' s as  a function  of  the  j ' ® 

®k  " ''k  ‘ic  ■ ‘Ic+l  = 2 ...  M (2.151) 

The  non-negativity  of  the  ® j ' s for  j =1  ...  M and  the  non- 
negativity of  oCj  for  j = 1 . . . M,  together  with  the  relation 
(2.149)  for  result  in  the  following  set  of  constraints  on 

the  waiting  times : 

0 < W2  £ 

0 < ^ ^ k = 2 ...  M-1  (2.152) 

“k41  - 0 

The  Jacobian  J of  the  inverse  transformation  between  waiting  times 
and  inter-arrival  times  described  by  (2.151)  is  the  following  matrix 


J = 


-1 

1 

0 


0 

-1 

1 


•-1/ 


(2.153) 


Since  J is  a lower  triangular  matrix,  the  diagonal  elements  are 
eigenvalues  and  hence  |det  J|  = 1.  Thus,  we  calculate 
p(Xj^  ...  •••  1^)  substituting  the  relations  for  0^ 

defined  by  (2.151)  into  (2.150)  and  combining  this  with  the  con- 
straints described  by  (2.152). 

r 

IT  B(x.)X  e 
• ^ ^ 

P(^l  • • • ^>^2  • ' • " fT 

M 

0 < < x^ 

0 < k = 2...M-l 


“^+1  “ 


(2.154) 


The  final  step  is  to  in 
service  requirement  x^. 


tegrate  out  the  dummy  variable  anO 


.M-1  . M-1  Ax. 

A B“(X)  , X-  1 

j: — ^ — ‘-  Tt  B(x.)e 

M i=l  ^ 


0 < CO2  £ Xi 


0 < I'ij+l  £ >^  + “V  k - 2 ...  M-1 


(2.155) 


B*(X)  denotes  the  single-sided  Laplace  transform  of  the  density 
B(x)  evaluated  at  X. 


We  observe  parenthetically  that  (2.155)  takes  a particularly 
simple  form  from  M/D/1  queue.  Since  the  service  requirements  are 

— 2 X 

deterministic,  B*(S)  = e .We  can  integrate 

over  all  the  impulses  to  obtain  the  joint  density  of  waiting  times 
conditioned  on  M.  Employing  the  explicit  formula  for  f^^,  the 
probability  of  M customers  being  served  in  a busy  period,  found 
by  solving  the  functional  equation  of  (2.139)  and  listed  earlier 
as  Eq.  (2.110),  we  find  that  p(c02  •••  is  given  by 


p(oJ2  ...  = M! 


0 < CO2  < X 

0 < k = 2...M-l  (2.156) 
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Evaluation  of  E(E  ,,  'M) 

n+k 


We  now  return  to  the  more  general  case  of  an  M/G/1  queue  and 
approach  the  calculation  of  the  mean  of  conditioned  on  M by 

deriving  the  joint  density  of  (x  ,o;  . . . cc  ) conditioned  on 

M.  We  use  this  joint  density  to  compute  the  distribution  func- 
tion of  Let  F ^ (t)  denote  the  probability  that 

^n+k IM 


C ‘ < T and  P 

n+k  — n 

n+k  (M 


(t)  be  the  probability  density  of  con- 


ditioned on  M.  Since  C is  defined  as  the  minimum  of  the  vari- 

n+k 

ables  •••  compute  the  probability  that 

< T as  one  minus  the  probability  of  the  complementary  event 
that  each  variable  is  greater  than  t. 


(t)  k Pr{c"  < tIM]  =1  -Pr[x„>T:oj  -T  >t:  ...  oj 


"n+k  \M 


(2.157) 


Now  we  calculate  the  density  of  by  differentiating  with 


respect  to  t. 


^n  " dr  ^(,n  , . 

n+k  \M  n+k|M^ 


(2.158) 


By  integrating  the  density  defined  in  (2.158)  we  can  calculate 

the  desired  conditional  mean  of  C^.,  . 

n+k 


■ V « •-  > s ♦ . 


' • *•  • V « Jt 


E(C“  ,„)  t C" 


■'n+k 


= n+k  \M 


J Tp  ^ (t)  dr 

T=0  ^n+k 


(2.159) 


We  obtain  the  joint  density  for  •••  ‘^n+k^  condi- 
tioned on  M by  integrating  first  over  the  variables  (0:2  •••  , 

<“n+k+l  “k>  <“1  •••  Vl-Vl  •••  V 

p(Xj^  ...  *'*  defined  by  (2.155).  For  n^l,  the 

suggested  integrations  are  performed  in  two  steps  and  the  results 

listed  below. 


P(^l  • • • • • • %-HclM^  J •*'  I 


X T + O!  , X + W . , 

n-1  n-1  n+k  n+k 


“2°*2  “■n‘*n 


CO 


h+k+l=0 


J 

p(Xj^  . . . ' ‘ ’ ’ ■ * *^^n+k+l  ’ * ’ 


n-1 

for  0<co,,  < S x.+x 
— n+1  — . , J n 
J=1 


0 < CO  , - < X . , + CO  , , 
— n+2  — n+1  n+1 


K = ^n+l 


- Z X ,o| 

j=k  ^ ) 


0<co.,  <x.,  ,+co,,  1 

— n+k  — n+k-1  n+k-1 


(2.160) 
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n 


Jt  • 


■nr  ' 


“n.k'"' 

n-1 


'n+1  ^n+1 


i-1 


1 


n.ax  -x^^qI 


^n+k-1  ^n+!'.-l 


X.  '^0  j = l ...n-1 
J - 


X , , = 0 

n+k 


’'m-1  ■ “ 


p (x,  ...  X,.  , a;  , , ...  a;  jj  ',M) 

^ ' 1 M-1  n+1  n+k 


dx. 


. dx  , dx  , , 
n-1  n+] 


. . dx 


M-1 

(2.161) 


A max 


Note  that  the  resulting  density  in  (2.161)  has  no  constraints 
remaining  on  any  of  the  variables  except  non-negativity.  The 
case  n=l  is  worth  distinguishing  since  in  the  final  density  one 
constraining  relation  remains  between  CO2  atid  x^. 

^l+k''‘‘*l+k  ^M-l"^^-! 

P(^X  ' ' * ^-l’^2  ■ ■ ■ ^ fk  ° I ' ■ ■ ' _ 

^'2+k=0  "1^=^ 

p(x^  ...  ...  c^’M)du^,  ...  du;2+k 

(2.162) 


96 


for  0 < 


0 £ <-1+^  £ ^ =-•  2 . . . k 


p(xj^,ui^  . . . |M)  J '''I 


^2~^2  ^1+k  ° ^M-1  ° 


p(x^  • • • ^m-1’‘^2  ■ ■ ■ ‘^+k  • • ■ *^^l+k  • • • ^^2 


(2.163) 


for  0 < cc.^  < x^ 


‘l  - "■“  l‘"'(,+l  - “t-Ol 


Hence,  the  distribution  function  F (t)  is  obtained  by 

c" 

applying  (2.157)  to  the  results  in  n-tk  |M  (2.161)  and  (2.163). 


n^l  1 - I I 


n n+1  n+k 


da;  ,,  . . . dx 
n+k  n 


■'n+k  IM 


1 - I I J 


Xi  = T iC^-T 


^+k* 


P ( X , o!2  • • • I M ) 


da;^^^  . . . dc02dxj^ 


(2.164) 


f 


p 


-,n 


(t)  could  now  L(j  obtained  by  di.f ferentiaLing  the  above 


n+k  ',M 

integrals  with  resp«..t  In  the  paraitiet(>r  t using  Leibni  tit's  rule. 

Simpler  results  for  possible  when  we  specialize 

the  preceding  calculat  Lons  to  the  case  of  an  M/D/]  queue.  Since 

the  service  requirements  are  all  the  same,  x.=x  and  C^,,  can  be 

1 n+k 

expressed  as 


^n+k  " ’ 


(2.165) 


where 


2 - '"i"  ...  . 


(2.166) 


By  inspection  we  write  the  density 
terms  of  the  distribution  function 
sity  of  z, 


of  C^,,  conditioned  on  M in 
n+k 

of  z,  F iw(t) , and  the  den- 
z I n 


p 

n+k ' M 


(t)  = 


(1 


^j^(x))6(t-x)  + U_^(x-t)P^(t) 

(2.167) 


There  is  a finite  probability  that  = x,  hence  the  impulse  in 

the  density  (2.167).  For  z < x,  = z and  hence  the  density 

of  P ^ (t)  , is  the  same  as  the  density  of  z for 

*"n+k  jM 
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0 < T < X. 


From  (2.167),  we  can  express  the  conditional  mean  of 


as 


VklH  ■ <1  - 


0 


(2.168) 


Hence,  for  an  M/D/L  queue,  the.  problem  of  computing 
iM  reduces  to  calculating  the  density  and  distribution  function 
of  z defined  in  (2.166),  conditioned  on  M.  We  need  the  joint  dis- 
tribution of  (u)  ,,  ...  o)  ,,  ) conditioned  on  M to  characterize  the 
^ n+1  n+k' 

distribution  function  of  z.  In  a manner  similar  to  the  deriva- 
tion of  (2.160)  from  (2.155),  we  integrate  out  the  unneeded 
from  the  density  p(o:2  •••  <^1M)  given  in  Eq . (2.156)  to  obtain 


P(i 


n+1 


‘^n+k  • 


M-1 

x+co„ 
X 2 

X+CO  , 

n-1 

'’<“,.+1  ■ 

■ • '«>  - <to> 

J f 

r 

J 

J 

il 

CM 

II 

'oT 

CO  =4> 

n ’^n 

‘‘n+k+l 

’■■'“k-i 

I ^“n+k+1 

^=0 


for  0 <:  CO  ,1  •-  nx 
— n+1  — 


0 < CO 

— n+6 


X + CO 


n+-t  - 1 


i = 2 . . . k 

^ max  foo^^j  - (n-t+l)x,0] 
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(2.169) 


I 


Hence,  by  the  s.ime  reasoning  tli.ii  led  to  (2.166), 
given  by 


- 1 - J 


X+U'  ,1  xbo!  ,,  1 

nx  n+1  n+k-1 


a;  ,-,-T  (Jj  ..,"T 

n+1  n+2 


’’^"n+l  "■  “n+k'”^ 


■^"n+k  • • • '^%+l 


Applying  Leibnitz's  rule  successively,  we  can  derive 


X+T 

’‘+“'n+2 

''+“n+k-l 

" f 

f 

•••  I 

^^■^’^^+2’  • • • ^n+k 

^n+2"'^ 

‘n+3"'^ 

“n+k"^ 

^‘^n+k  •'*  ^^^^+2 

nx 

X+T 

’‘■"“n+k-1 

+ J 

J 

•••  f 

P(%+l'"’"^n+3  •••  %+k 

to  ,,  =T 
n+1 

‘n+3"'^ 

Si+k=T 

^‘^n+k  • • • ‘^'*n+2 

x4co 

nx  n+1 


+L  .i=T  to  .o  = T tO_.,-T 


"n+1  ■ n+2 


p(u,;^+l  . . . +Vi+k-2’'^’^n+k'^^ 


d CO  I , ...  d CO  , 1 

n+k  n+1 


x+u:  , , 
nx  n+l 


n+1  ■ “^n+2 


n+k- 2 


+L  I 1 T t0_  I r,  T t0„  .1  T 


n+k-1 


pCoj^^j  ...  \M) 


n+k-1 


(2.171) 


In  the  next  few  paragraphs  we  show  liow  to  compute  explicitly 


C , We  were  not  able,  however,  to  obtain  explicit  terms  for  a 

n+1 


general  C and  therefore  have  to  rely  on  (2.173)  when  other 


m 


C”lM  will  be  needed.  For  IM,  z becomes  the  single  variable 

m ‘ n "tI  ' 

cc  Note  that  to  compute  the  conditional  mean  in  (2.168)  we 

n+1 

only  need  ^ ^ ^ Hence,  the  approach  we  take  is 


to  derive  the  explicit  form  for  P 


Ui 


n+llM 


(t)  on  0 < t < X.  From 


(2.169),  we  write  the  density  for  P 


CO 


(t)  on  0<T<x  as 


h+llM 


^n+1  IM 


(t)  = 


M! 


(Mx) 


M-f  -■ 


X -t  u\,  X + a;  1 x+r 

2 n- 1 


u.'2~0 


'"n=0 


n+2 


■I 

■®  “w-o 


' du:  . odio 
n+2  n 


• OCl 


n-1 

Integrations 


M-n-1 

Integrations 


(2.172) 

Note  that  the  above  integral  consists  of  two  sets  of  inte- 
grations defined  by  the  brackets.  We  define  the  following  set  of 
iterated  integrals 


X+CO 

I^^^(a')  = J I^(o)' )du>'  with  Ij^(co)=l 


(2.173) 


CO 


j-0 


Hence,  we  can  express  I’ 

I.  's. 

J 


CO 


(t)  for  0 • T < X in  terms  of  these 


hillM 
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(2.174) 


Hi 


1 


M-  1. 


P (t)  = M!  (^-)  1 (0)1„  (t) 


0 < T < X 


Examining  a few  of  the  I.'s,  we  assune  that  I (u))  can  be 

J n ' ^ 

represented  as  a power  series  in  (x  + u,')  with  powers  up  to  (n-l)-s 
order . 


n-1 


I (co)  = S 


n 


i=0 


(2.175) 


,(1)  = 


To  be  consistent  with  I^(co)  = 1,  = 1.  Employing  (2.173)  we 

can  derive  a set 


's. 


ference  equations  relating  to  the 


n-1 


n("+i)  - E i ( i 

q q q-o 


q=l  ...  n 


(2.176) 


It  is  verified  in  Appendix  A that  the  solution  for  the  of"^'s 
are  given  by 
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\ .7  V - 


1 < 


I 


for  n > I n 


(n)  _ 


n-2-i 

X ill-Jl 

(n-1-1  '!(i-1)! 


i = 1 ...  n-1 


i = 0 


for  n = 1 = 1 


(2.177) 


Employing  (2.168),  the  expression  for  the  density  P 


n+1  IM 


given  in  (2.174),  and  the  definition  of  the  ^ j ' s in  (2.173), 
|M  may  be  expressed  as 


M-U 


VilM  - (1  - M!(^)  I„(0)I„.„^l(0)). 


M-1- 


(2.178) 


By  looking  at  P (t)  for  n = 1 we  can  obtain  an  expression  for 

‘Si+l  |M 

Ij^(O)  since  the  density  for  a;^  must  integrate  to  1 over  the  interval 
[0,x]  . Hence,  Ij^^(O)  must  be  given  by 


I (0)  = — 

^ M! 


(2.179) 


To  specify  the  only  remaining  unknown  in  (2.178),  we  define 


m ^ ''  tI  i't)  dr 

n - n 


(2.180) 


• V 


We  evaluate  this  integral  by  substituting  our  explicit  form  foi 


I (t)  defined  by  (2.175)  and  (2.177) 
n 


m = 
n 


n = 1 


(2.181) 


V 

Putting  together 
expressed  as 


n-1 

L 


n+1.  ,.n-2-i  i , . 

X (n-1) y _X- 

(n-l-i)!(i-l)!  t+2 


n 1 


(2.179),  (2.181),  and  (2.178),  may 


be 


(1  - Mf( 


M-1 


n!  (M-n+1)!  ' 


M-1 


(nx) 


n-1 


n! 


“Vi- 


n 


(2.182) 


Bias  Calculation  for  an  M/P/l  Queue 

Having  defined  the  calculation  of  \M=i  and  hence  Q(i)  , 

we  can  investigate  the  behavior  of  the  asymptotic  bias  b expressed 
as  a power  series  in  p by  Eq.  (2.144).  For  simplicity,  we  start 
with  an  M/D/1  queue.  - 0 in  (2.136),  since  the  variance  of  the 

service  time  density  is  zero.  Since  each  ^ function  of 

(x,o)^_^^  ...  summation  of  for  a given  busy 

period  with  M customers  can  be  expressed  as  different  functions 
of  (x,t«:^  ...  c»i^)  over  regions  in  (u>2  •••  space.  Hence  Q(M) 
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can  be  expressed  as  the  sum  of  a set  of  integrals  of  functions 

of  (x,W2  •••  over  regions  in  (u:^  ...  uj^)  space  weighted  by 

1 M-1 

the  density  p(u;2  •••  ^ specified  in  Eq.  (2.156). 

Since  the  density  p(ti,'2  •••  is  a constant,  independent  of  p, 

Q(M)  will  also  be  a constant  with  no  dependence  on  p.  f^,  the 
probability  of  n customers  being  served  in  a busy  period,  is 

j-  p^  Hence,  expanding  f^  as  a power  series  in  p sub- 

stituting the  result  into  the  expression  for  p in  (2.142),  and 
changing  the  order  of  summation,  we  can  derive  an  explicit  formula 
for  a^,  the  coefficient  of  in  a power  series  expansion  of  p 

in  powers  of  p.  Therefore,  on  the  basis  of  (2.143),  to  show  tha*" 
the  estimation  procedure  is  unbiased  up  to  the  -t“th  power  in  p, 
we  have  to  show  that 


ot_ 


n 


n+2 

r 

i«=2 


i! (n+2-i)! 


Q(l) 


(n+2)(n+lj 


n 


t-1 


(2.183) 


By  a change  of  variable  n * n’-2  and  some  manipulation,  verifying 
conditions  (2.183)  may  be  reformulated  as  checking  the  following 
recursion  for  the  Q(i)*s: 
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Q(n')  = 


!lJ (Hl(n'-l)  ) 


n! 


(n')"'-l 


n'-i.n'-l 


n'-l  , . 

s iil 

i=2 


i! (n'-l) ! 


Q(i) 


n ' =3  ...  -C+1 


with  Q(2)  “ ^ X 


(2.184) 


Employing  the  preceding  results  for  an  M/d/1  queue,  we  pro- 
ceed to  show  that  the  asymptotic  bias  b only  contains  powers  of  p 
greater  than  third  order.  For  -1=3,  relation  (2.184)  dictates 
that  Q(2),  Q(3) , and  Q(4)  assume  the  following  values; 


Q(2)  = J X 

Q(3)  = I X (2.185) 

Q(4)  - X 


By  definition  of  the  Q(i)'s  given  in  (2.134),  Q(2),  Q(3),  and 
Q(4)  are  given  by 


We  can  compute  all  C^|M's  of  the  form  \M  by  employing  the 

general  formula  derived  as  Eq.  (2.182).  To  calculate  other 
C^\M's  we  appeal  to  the  procedure  outlined  in  Eqs.  (2.165)  through 
(2.171).  The  results  of  the  calculations  implied  in  (2.186)  are 
listed  below. 


7^2  ,,  _ 19 

(2.187) 

By  using  the  results  of  (2.187)  in  (2.186)  it  is  easy  to  see  that 
the  numbers  for  Q(2) , Q(3),  and  Q(4)  are  consistent  with  the 
values  listed  in  (2.185). 

For  the  sake  of  clarity,  we  go  through  the  derivation  of 

"”2 

C^|4  as  a sample  calculation.  The  remaining  C^lM's,  those  not 

2 

covered  by  relation  (2.182),  are  calculated  in  Appendix  B. 
is  defined  as  follows: 


cl\2  = i X 

C2I3  = I X 

C2l3=f  X 
C2'14=|Jx 

C3I4  = I X 

C4I4  =f|x 
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with 


A - ■ -s 

- nun  lx.zj  , 


(2.188) 


z = min 


(2.189) 


the  joint  distr^.bution  of  ( 0)2 > » co^ ) conditioned  on  M=4  is  given 
by  Eq.  (2.156)  as 


P(a!2,a.-3,a;^  1M=4)  = 3 

8 X 


(2.190) 


0 < 0^2  < X 


0 < 0C3  < X + 0^2 


0 < oj,  < X + uJo 

- 4 - 3 


Hence,  the  joint  distribution  of  (c02,to^)  conditioned  on  M=4  is 


given  by  Eq.  (2.169)  as 


p (0:3,0:^  |M=4)  = ^ J 


dco,  . (2.191) 


OJ. 


>2  = max  [0,a;3-x] 


0 < CO3  < 2x 


0 < 


- 4 < X + 0^ 


More  explicitly,  (2.191)  means  that 
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•»r  «:t  .>  V - 


3 

pCoJ^.O!^  lM-4)  = 2 

8x 

and 

p(tcv  ,UJ,  \M-4)  = 3 (2x  - 

^ ^ 8x^ 

u>3 ) • 

0 < < X 

X < 0)3  £ 2x 

0 < < X + 

0 < tO/  < X + to- 

— 4 — j 

(2.192) 

Hence,  by  (2.170)  the 

distribution  function  of  z conditioned  on 

M=4  and  for  z in  the 

interval  [0,x]  is  given  by 

X 

X + OJ^ 

1 - I 

r 9 dto,  dw. 

Wo  = T 

0 < T < X 

C0^=T 

o X + 

2x  3 - 

- J f 

0X^=X  C0^=T 

J 

8x^ 

(2x  - 0O3)  800^8003  . 

(2.193) 

Differentiating  (2.193)  with 

respect  to  t and  doing  the 

integra- 

tions  yields 

15  3 

^ ■ 8x2 

(2.194) 

0 < T < X 

From  Eq.  (2.168),  we 

compute 

—2 

C^\4  as  follows: 

X X 

C4IA  = x(l  - f P2|m=4^^>  I ’■^z\M=4^''^  dT  = ^ X 

• 

0 

0 

(2.195) 
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0t  • ♦ 

•%  4-- 

Bias  Calculation  for  an  M/M/1  Queue 


We  now  examine  the  behavior  of  the  asymptotic  bias  b for  the 
case  of  an  M/M/1  queue.  The  service  requirements  are  now  expon- 
ential random  variables  with  parameter  ji.  Hence,  B(x)  and  B*(s) 
are  given  by 


B(x)  = 


3*(s) 


(2.196) 


s + p 


Since  the  mean  service  time  is  — and  the  variance  of  the  service 
1 2 

requirement  is  — the  defined  in  Eq.  (2.136)  is  one.  The 
functional  equation  of  (2.139)  may  be  solved  to  obtain  F(z) 
and  hence  f^,  the  probability  of  i customers  being  served  in  the 
busy  period 


£ . 1 (21-2,  y... 


(2.197) 


To  characterize  the  asymptotic  bias  b we  need  Uy  the  coef- 
ficient of  in  a power  series  expansion  of  p,  which  is 

defined  in  (2.132)  as  a summation  of  the  products  Q(i)f^  for 
i = 2 ...  «>.  We  contend  that  Q(i)f^  can  be  expressed  as  a power 
series  in  p whose  terms  are  at  least  (i-l)-st  order  in  p.  Hence, 

to  find  a.  we  need  to  collect  coefficients  of  p^^^  in 
J 

Q(2)f2  ...  Q(j+2)fj^2‘  reason  for  Q(i)f^  being  expressible 
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W ■ ■'  •V 


as  a power  series  in  p with  terms  of  at  least  (i-l)-st  orders 


follows  from  Eq.  (2.155)  for  p , CO2  • • • \M)  . Q(M) 

represents  the  mean  of  a random  variable  which  is  expressed  as  a 

summation  of  C^'s,  each  of  which  is  a function  of  (x^^  ... 

o:^  . . . <^)  . Hence,  Q(M)  could  be  expressed  as  an  integral  of 

various  functions  of  (xj^  ...  weighted  by 

p(x^  ...  •••  over  (xj^  ...  space. 

From  the  density  given  in  (2.155),  this  integral  for  Q(M)  will 

M-l 

have  a leading  factor  of  X / Hence,  Q(M)fj^  can  be  expressed 

M-1 

as  the  product  of  p and  an  appropriate  integral.  This  shows 
that  Q(M)fj^  can  be  represented  as  a sum  of  terms  which  are 
(M-l)-st  order  or  higher  in  p.. 

We  now  proceed  to  compute  Q(2)  and  Q(3)  so  we  can  calculate 
Oq  and  and  then  the  coefficients  of  the  first  and  second  power 
of  p in  the  asymptotic  bias  b defined  in  Eq.  (2.144).  Q(2)  is 

equal  to  \2  and  is  given  by 


C2  “ min  {Xj^,c02}  . 


(2.198) 


From  (2.155),  the  joint  density  of  (xj^,a'2)  conditioned  on  M=2  is 
as  follows: 


P(x^,o!2  1M=2)  = — B(Xj^)e 


(2.199) 


0 < o:^  < Xj^ 


S 


Hence,  by  (2.164)  the  d Ls  trilmt  i (iii  function  of  conditioned  on 
M=2  is 

‘ ■ I R(x^)e  . (2.200) 

:<1-T  lc^-t 

Differentiating  with  respect  to  ^ by  applying  Leibnitz's  rule, 
we  obtain 


V,|M.2<->  <iXl  (2-201) 

Xl  = T 

Employing  (2.201),  we  can  express  the  desired  conditional  mean  as 
a double  integral.  Changing  the  order  of  integration,  we  obtain 
.he  following  result  for  an  M/G/1  queue: 

CD  “X  X 

Q(2)  = - I j ^2^  ^B(x^)dx^  (2.202) 

x-j^=0 

For  an  M/M/1  queue,  (2.202)  specializes  to 


0(2)  - O2U  " \ i 


4 f 

(1+p)  2 


(2.203) 


”1—1  —2 

Ne.xt,  we  calculate  <)l3)  as  the  sum  of  C,,|3,  C^'3,  and  C^13. 
iiie  joint  di-nsity  of  (Xj^.^'^)  conditioned  on  M=^3  is  .specified  by 
Eqs.  (2.155;,  (2.162),  and  (2.163)  as 
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X , + 


I I 


-A.  (x^+x^) 

Bf  B(x2>e  doo^dx^  . 


0 < 0^2  < Xj^ 


x^  =0  u.',=0 

I ) 


(2.204) 


Performing  the  inner  intc^gratLon  over  u:^  and  changing  the  order 
of  the  integration  over  ihe  v.iriable«  ^2  and  Xj^  specified  in 
(2.164),  the  distribution  function  of  conditioned  on  M=3  is 
given  by 


CO  9 


1 ,M„3(0  - 1 - — J J f (x2+a-2)B(xpB(x2> 


-A (x^+x^) 


qiM==3 


tu^^T  X2=0 


dx^dxj^do;^  . 


(2.205) 


Differentiating  (2.205)  uith  J.iebnitz's  rule  we  obtain  ^j,j,,-^(t) 

Based  on  this  P^iiw  -iCt)  we  can  write  the  following  integral 

C2 1M=3 

expression  for  the  conditional  mean: 


2 CO  OD  00  ■■  A (Xj^+X2) 

cha  = ^ J T J J (x2+r)B(xj^)B(x2)e  dx2dx^dr 


T=0  Xj^  = 7 X2”0 


(2.206) 


For  ar  M/M/'l  queue  wo  employ  the  B(x)  specified  by  Eq . (2.196) 
and  obtain  the  result 


7:1,  ^ . 1 P _ 1 

(.„13  = 3 - — 7 7-  . 

^ ^ (1  t-p)^  ^3 


(2.207) 


W • -'V 


By  Eqs.  (i.155),  (2.160),  anc  (2.161),  the  joint  density  of 
(x2,o!j)  ccnditiimed  on  M = 3 is 

P(X2,C03  1M=3)  = J J 

=max  (a)3-X2,0]  003  = max  [u;3~X2,01 

-X(x,+x») 

B(x3)B(x2)e  ^ dw2dx3  . (2.209) 


1^ 

■< 

■i 


▼ 


Doing  the  inner  integration  over  (^2*  and  noting  the  form  for 

F 2 dictated  by  (2.16A),  differentiation  of  F „ (t) 

''3\M=3  €3111=3 

w.th  r?spect  to  r by  Leibnitz's  rule  will  yield  the  following 

integral 


CO 


p 2 _(r)  - I 

C3lM=J  X2 


P(x2,TlM=3)dx2  + f P(T,ui3  lM=3)du.'3 


(2.210) 


Employing  (2.209)  and  (2.210),  the  desired  conditional  mean  may 
be  expressei  as  the  following  integral : 


C?\3 


« -AX  -Xx 

r J E(<2>e  " x^»(Xj^)e  ^ dxj^dx^dr 


T“0 


J TB(i)e^^ 

T=0 


ao  oc* 


Xj-O 


-Xx, 


J [ Xj^  - (co^-t)  ] dx^dto^dr 


<0^=T  Xj^  = a'j-T 


(2.211) 


Using  the  8(x)  corre Jponding  to  an  M/M/1  queue,  we  obtain 


= 2 - ^ 
3 ' 


7-1. 


(2.212) 


Our  final  calculation  is  of  C^p.  The  quantity  is 
definiid  by 


“ min 


(2.213) 


The  joint  density  of  (Xj.a)2,u.j)  conditioned  on  M=3  is  determined 
by  Equ.  (2.155),  (2.-62),  and  (2.163). 


P(Xj^, (02*^3  \M=3)  * J 

0 < oj.,  < Xj^ 


3(Xj^)  B(x2)e 


-X  (x^-h<2) 


dx. 


<2  ' niax  (.uJ2-u!2 ,0  ; 


(2.214) 


Employing  (2.164)  to  derive  tne  distribution  function  \M=3^^^ 
and  breaking  up  the  Integrati  m over  to  get  rid  of  the  "max" 
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in  the  limits  of  the  integration  over  we  obtain 


X OJ 

1 2 


I J J I B(x^)B(x2>e 

Xj^=T  co2=t  oj3=t  X2=0  dx2 Jw3du;2dx^ 


-X (X3+X2) 


A 1 

CD  CO  OD 

r r •'  r 

J J J J 


B(xj^)  B(x2)e 


-X(x3+X2) 


X,=T  OJ„  = T 0)o  = U;„  Xo  = Wo-a!„  J , J . 

1 2 3 2 2 3 2 dx2dc03du)2dx3  . 


(2.215) 


Differentiating  (2.215)  by  applying  Liebnitz's  rule  we  can  derive 

P^liw  o(t)  and  from  that  a single  integral  expression  for  C^\3. 
C3iM=3  -J 


-1.  Q 

C3|3  = 


CD  CO 


-X (X3+X2) 


J T J J J B(x^)B(x2)e  dx2dG02dx3dT 


T=0  X3=T  0)2=T  X2-O 


2 CD  CO  CD  00 

+ L-BiiXi  J.  ^ J J.  J. 

3 

T=0  x^=T  u:3=t  X2=oo3~t 


B(x3)B(x2)e 


-X(x3+X2) 


dx2du,'3dx^dT  (2.216) 


Using  the  B(x)  for  an  M/m/1  queue  given  by  (2.196),  we  find 


c\\3  = 5 - j-  . 

^ ^ (1+p)^  ^3 


(2.217) 


Now  we  can  calculate  the  coefficients  in  the  asymptotic  bias 
b for  the  first  and  second  order  terms  in  p.  Employing  (2.203), 
(2.207),  (2.212),  and  (2.217),  we  compute  Q(2)  and  Q(3)  as 


•4  - 


Q(2)  = 


1 1 D _1_ 

2 u . sA  f. 


^ (1+p)  2 


(2.218) 


Q(3)  = 


10  - 


^ _L 

^ (1+p)^  ^3 


ttQ  is  defined  by  the  coefficient  of  p in  Q(2)f2*  is  deter- 


mined by  collecting  powers  of  p in  Q(2)f2  and  Q(3)f2*  Hence, 
we  obtain 


1 1 
2 


(2.219) 


a,  = 


8 


Therefore,  from  (2.144)  we  can  represent  the  asymptotic  bias  b as 


, 11  . 11  1 2 . 3. 

■ 2 * T ^ + 0(p  ) 


(2.220) 


O(p^)  denotes  terms  that  are  third  order  or  higher  in  p.  This 
shows  that  the  algorithm  is  biased  for  M/M/1  and  contains  terms 
of  all  orders  in  p,  in  contradistinction  with  the  M/D/1  case 
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where  we  showed  that  the  algorithm  is  asymptotically  unbiased  U(> 
to  third  order  in  p.  In  fact,  we  are  confident  that  for  an 
M/d/1  queue  the  algorithm  is  completely  unbiased,  although  we 
were  not  able  to  show  this. 


2 . 7 Cramer-Rao  Bound  for  Customer-Removal  Algorithm  in  Case 

of  M/D/1  Queue 

Now  we  derive  a Cramer-Rao  bound  for  the  variance  of  any 
unbiased  estimator  of  the  delay  gradient  that  works  with  the 
same  observations  as  the  customer-removal  algorithm  in  the  case 
of  an  m/d/1  queue.  Asstuning  the  service  requirement  is  a known 
variable  x,  the  observations  which  the  customer  removal  algorithm 
processes  to  derive  the  c'^'s  and  form  the  estimator  are  the  M.'s, 
the  number  of  customers  served  in  the  i-th  busy  period,  and 


(w, 


(i) 


Xi) 


'2  * * • ‘‘m  ^ ^ waiting  times  of  customers  in  the  i-th  busy 

i 

period.  Hence,  we  define  our  observation  vector  Y'  as 


Y'  A 


(1) 


(N) 


(N). 

N 


(2.221) 


Since  waiting  times  and  the  number  of  customers  served  per  busy 
period  are  statistically  independent  from  one  busy  period  to 
another,  the  joint  density  of  Y'  may  be  expressed  as 


P(Y') 


IT  P(M 
i=l 


1 
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(2.222) 


-•> 


We  compute  the  joint  density  of  ...  by  breaking 

it  into  the  product  of  ...  ) and  Pr  (M.).  Hence, 

employing  Eq.  (2.156)  for  the  joint  density  of  waiting  times 
conditioned  on  the  number  of  customers  served  in  a busy  period, 
expressing  the  densities  in  terms  of  the  parameter  p = Xx,  and 
applying  (2.222)  we  obtain 


P(Y’  Ip) 


/ N ' 

S M.  - N 
p\i=l  ^ ! 

r \ 

L M.  > N 

x\i=l  ^ / 


I ^ 

L 

\i=l 


0 < < X i = 1 . . . N 

0 £ ^ k = 2 ...  M^-1 


(2.223) 


Lei.  ting  y denote  the  delay  gradient,  y and  p may  be  related 

by  (2.112).  Employing  (2.113),  we  can  derive  — Using 

by 

the  expectation  in  (2.114)  and  the  general  formulation  given  in 
(2.102)  for  a Cramer-Rao  bound  on  the  variance  of  any  unbiased 
estimator  of  a parameter  4 based  on  an  observation  vector  R,  we 
find 

Var  (y  - y)  > x^'  ^ r • (2.224) 

N(1  -p)^ 
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We  can  shed  light  on  the  reason  for  the  p in  the  numerator 

of  the  bound  for  the  customer-removal  algorithm  as  compared  to  the 
2 

p in  the  numerator  of  the  bound  for  the  customer-addition  algor- 
ithm by  examinixig  the  maximum  likelihood  estimation  procedures 
that  follow  from  the  observation  vectors  Y and  Y'  employed  by  each 


p 


N 

L 

i=l 


N 

L M.  - 1 

i=l  ^ 

N-1 

M.  + L I.  lx 


(2.228) 


The  estimator  that  follows  from  the  observation  vector  Y'  employed 
by  the  customer-removal  algorithm  is  given  by  (2.226)  as 


^ = 1 - . (2.229) 

S M 
i=l  ^ 

Each  estimator  for  p can  be  shown  to  be  asymptotically  unbiased 
and  together  with  (2.227)  must  determine  an  asymptotically 
unbiased  estimator  of  the  delay  gradient. 

Since  the  delay  gradient  estimator  specified  by  using  (2.228) 
in  (2.227)  uses  idle  period  information  in  addition  to  the  number 
of  customers  served  in  each  busy  period,  we  might  expect  that 
asymptotically,  the  variance  of  this  estimator  will  be  smaller 
than  that  for  the  estmmator  that  follows  from  using  (2.229)  in 
(2.228),  Hence,  it  is  not  surprising  that  the  Cramer-Rao  bound 
for  the  customer-addition  algorithm  is  always  smaller  than  che 
corresponding  bound  for  the  customer-removal  algorithm. 
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2 . 8 Derivation  and  Realization  in  Flow  Diagrain  Form  of 
Time-Contraction  Alj^orithm 


In  our  third  algorithm,  we  simulate  an  increase  in  rate  6x 
by  a linear  contraction  in  time  scale.  Assume  that  i denotes 
the  time  of  arrival  of  the  j-th  customer  in  the  i-th  busy  period 
relative  to  the  beginning  of  the  observation  inteirval.  We 
define  a new  set  of  shifted  arrival  times  by 


^*(i)  = (1  _ • (2.230) 

Since  6X  represents  an  infinitesimal  change  in  rate,  we  can  choose 
it  sufficiently  small  so  none  of  the  busy  periods  are  shifted 
onto  other  busy  periods  by  the  time  contraction.  A simple  suf- 
ficient condition  on  6X  so  that  the  redefinition  of  arrival  times 
in  (2.230)  produces  no  "interactions"  between  busy  periods  is 
given  by 


X 


< min 

E i 


(2.231) 


T refers  to  the  observation  period  length  and  I.  denotes  the 
E J 

duration  of  the  j-th  idle  period. 


We  now  consider  the  increment  in  system  time  that  comes  from 

each  customer  arriving  a little  earlier  and  hence,  waiting  a 

little  longer.  The  waiting  time  of  the  j-th  customer  in  the  i-th 
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busy  period  is  defined  by 


. 'z  - (T<i>  - r<^>)  . 

J <,-1  <■  J ^ 


(2.232) 


If  we  substitute  the  shifted  arrival  times  given  by  (2.230)  into 
(2.232),  we  can  relate  the  new  waiting  times  to  by 


*(i)  = (i)  + ^ rT-(i)  . /i)i 

coj  - ojj  + ^ j . 


(2.233) 


We  define  a new  variable  Tj 


R(i) 


as  the  arrival  time  of  the  j-th 
customer  in  the  i-th  busy  period  relative  to  the  start  of  that 
busy  period.  Hence,  the  additional  system  time  over  N busy  per- 
iods that  follows  from  Eq.  (2.232)  can  be  expressed  as 


■<  ^ 


where 


f h 

^ i=l  ^ 


Z.  L 

X 


M. 

1 

L r 
n=2 


R(i) 


n 


(2.234) 


M.  = 1 

1 


(2.235) 


1 


A second  contribution  to  the  increment  in  system  time  follows 
from  the  fact  that  our  Lime  contraction  procedure  shifts  the 
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right  edge  of  the  observation  interval,  leaving  a gap  of  duration 
— in  which  additional  customers  could  arrive.  Since  an 
average  of  XD^  total  delay  is  accumulated  per  unit  time,  wheie 
is  the  average  total  delay/customer,  the  average  increment  in 
system  time  is  given  by 


(2.236) 


Employing  (2.234)  and  (2.236),  we  obtain  the  following 
expected  increment  in  system  time  conditioned  on  the  queueing 
record. 


E(6S (Queueing  Record)  = 6XT  D + L Z.  (2.237) 

^ C A • 1 X 

1=1 


Since  both  X and  D^"  are  unknowns,  we  use  the  fact  that  asymp- 
totically, by  the  law  of  large  numbers 


L M. 
i=l  ^ 


N "i  N 

L L Si^V  E M. 


i=l  j=l 


i=l  ^ 


(2.238) 


(2.239) 


Sj  denotes  the  system  time  of  the  j-th  customer  in  the  i-th 
busy  period.  Substituting  relations  (2.238)  and  (2.239)  into 


1 


r 


k 

( 


(2,237)  and  normalizing  by  6XT,,  we  obtain  the  following  delay 

h 

gradient  estimator; 


N ^i  ...  N 

L L sr''  Z Z. 

D'  = ^ ^ (2.240) 

L M.  L M 

i=l  ^ i=l  ^ 


Comparing  Eq.  (2.128)  and  Eq . (2.240),  we  see  that  the 
customer-removal  and  time-contraction  estimators  have  an  iden- 
tical first  term  corresponding  to  an  estimate  of  D^.  They  differ 
only  in  the  second  term  where  the  customer-removal  algorithm 
employs  p^^  defined  by  (2.129)  and  the  time  contraction  algorithm 
uses  specified  by  (2.235).  Hence,  flow  diagrams  for  both 
algorithms  are  very  similar.  The  variables  that  appear  in  the 
flow  chart  are  the  same  as  those  for  the  customer-removal  algor- 
ithm defined  in  Fig.  2.5  except  that  we  no  longer  need  the  ' s 
and  S is  redefined  as  the  sum  of  all  the  system  times  and  rela- 
tive arrival  times  of  customers  in  the  current  busy  period. 
Dropping  the  indices  for  the  busy  period,  S is  given  by 


S 


L (x.  + CO.)  + E . 

j=l  J J j=2  J 


(2.241) 


Noting  that  cOj=0  and  expressing  co^  as  the  sum  of  the  j-1  preced- 

R 

ing  service  times  minus  the  arrival  time  Tj  of  the  j -th  customer 


W ' * V 
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relative  to  the  start  of  the  busy  period,  we  obtain  the  following 
expression  for  S since  the  contribution  due  to  the  relative 
arrival  times  is  canelled: 

M M i-1 

S = L X.  + S L X . (2.242) 

i=l  ^ i=2  1=1 

We  use  (2.242)  in  the  flow  diagram  for  the  time-contraction 
algorithm  presented  in  Fig.  2.7. 

The  expression  for  S given  in  (2.242)  suggests  an  alterna- 
tive form  for  the  estimator  in  (2.240).  Changing  the  order  of 
summation  in  the  double  sum  in  (2.242)  and  realizing  that  S is 
the  contribution  to  the  numerator  of  (2.240)  corresponding  to  a 
given  busy  period,  we  obtain 


where 


B. 

J 


Z 

j=l 


M. 

J 


M. 

B.  A L (M.  - t + l)xjj^ 
J ~ t=l  J ^ 


(2.243) 


(2.244) 
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Figurr  2.7  Flow  Diagram  for  Time-Contraction  Algorithm 
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1 1.>  ..  X*. 


r ' 1 

2 . 9 Asymptotic  Bias  Calculation  for  Time-Contraction  AlK^rithm 
• in  the  Case  of  M/g/I  QueucT 

To  investigate  the  asymptotic  properties  of  the  time- 
contraction  algorithm  we  can  employ  the  form  for  the  estimator 
given  in  Eq.  (2.240)  or  (2.243).  We  first  examine  Eq.  (2.243), 
since  the  asymptotic  unbiasedness  of  the  algorithm  for  an  M/D/1 
queue  follows  readily.  Exchanging  the  limit  and  expectation 
operations  and  appealing  to  the  law  of  large  numbers,  we  obtain 


lim  ED'  = E lim  D'  = ® . 
00  ^ 00  ^ 


(2.245) 


For  an  M/D/1  queue  x^"^  ^ = x and  we  evaluate  B by  conditioning  on 
M=j  and  then  taking  the  expectation  over  M.  Hence,  B is  given  by 


■ 


B = Ej^E(B|M=j)  = j X + M . (2.246) 


By  Eq.  (2.245),  the  mean  of  the  estimator  asymptotically 
approaches 


y x(  ^ + 1 1 . (2.247) 

' M 


Employing  M/D/1  formulas  for  M and  M given  in  Eqs.  (2.98)  and 
(2.99),  we  obtain  the  final  result  that 
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1 * * 


V?  t;> 


(2.248) 


lim  ED  = X + 
N “ 


_£2L 


-fi2L 


2a-p)' 


The  expression  in  (2.248)  is  identical  to  that  in  (2.33)  for  the 
delay  gradient  of  an  M/D/1  queue. 

We  can  extend  the  argument  of  the  preceding  paragraph  to  an 
M/G/1  queue  for  calculation  of  the  asymptotic  form  of  the  esti- 
mator as  a poiver  series  in  p . To  evaluate  the  expectation  of  B 
conditioned  on  M we  need  p(x^  ...  IM)  since 


E(BlM)  = E(Mx^  + (M-1)x2  + ... 


(2.249) 


By  Bayes ' rule 


p(Xj^  ... 


p(M\Xj^  ...  Xj^)p(Xj^  ...  Xj_j) 


M 


(2.250) 


We 


calculate  p(M\xj^  ...  x^^)  as  the  probability  of  the  region  in 
(0j^  ...  6j^)  space  that  corresponds  to  M customers  with  service 
requirements  (x^^  ...  Xj^)  being  in  the  busy  period.  By  the  con- 
straints on  the  0 j ' s specified  in  (2.148)  and  (2.149)  for  M cus- 
tomers to  be  served  in  a busy  period  and  since  the  unconditional 
joint  distribution  of  the  Oj's  is  a product  of  M exponential  den- 
sities with  parameter  X,  p(M(xj^  ...  x^^)  is  given  by 
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p(Mlx^  . . . x^)  = J 


e^=o  02=0 


•-1  1 ^ 
1=1  1=1 


0^=0 


L X. 

i=l  ^ 


E 0. 

i=l 


X^e  ^ ^ d0M  • • • 

M i 


(2.251) 


Doing  the  final  integration  over  0^^  and  combining  with  (2.250) 
we  obtain 


M 1=1 


(2.252) 


where 


g(Xj^  . . . 


^1  x^+X2"®i^ 

I I 

0^=0  02=0 


E X. 
i=l  ^ 


Vi=° 


M-2 
E 0. 


^®M-1  • • • 


M 1 


(2.253) 


Hence,  we  can  calculate  B as 


B = E (B\M)f 
M=1 


M 


(2.254) 


From  the  (2.252),  we  can  see  that  (BlM)fj^  will  be 

representable  as  a power  series  in  p with  terms  of  at  least 
(M-l)-st  order.  Forming  B(M  by  multiplying  B with  1-p,  we  can 
collect  the  powers  of  p in  the  asymptotic  form  of  the  estimator. 

We  have  demonstrated  that  contrary  to  what  happens  for  an 
M/D/1  queue,  for  the  more  general  case  of  an  M/G/1  queue,  analysis 
of  the  expression  given  by  (2.243)  for  the  time-contraction  esti- 
mator does  not  yield  an  explicit  closed  form  expression  for  the 
asymptotic  mean  of  the  estimator.  Hence,  rather  than  pursuing 
the  calculation  outlined  by  Eqs . (2.245)  through  (2.254),  we 
examine  the  form  of  the  estimator  in  (2.240).  Comparing  this  with 
(2.128),  we  note  that  plays  the  same  role  as  p^  in  the  customer- 
removal  algorithm.  Therefore,  the  asymptotic  bias  of  the  time- 
contraction  algorithm  for  an  M/G/1  queue  is  specified  by  Eq. 

(2.134)  where  is  defined  as  the  coefficient  of  in  a 

power  series  expansion  of  Z.  Z is  expressed  using  (2.235)  as 


E E 
M=2 


^ R 
E r^lM 

i=2  ^ 


“■M 


(2.255) 
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In  the  next  few  paragraphs  we  develop  results  useful  in  cmI- 
/ M O \ R 

culating  E L t . iM  . Since  t.  is  equal  to  the  sum  of  the  first 
\i=2  ^ ^ \ 


(i-1)  inter-arrival  times,  the  problem  of  computing  E L t.\M 

\ i=2  " 

can  be  restated  as 


L T . \M  = E L e . (M- j ) \M  . 
i=2  ^ j=l  J 


(2.256) 


To  evaluate  the  second  expectation  we  need  p(6|^  ...  IM)  . 

Integrating  out  0^^  and  from  p(x^  ...  ■*'  ®M  defined 

in  Eq.  (2.150),  we  obtain 


P(xi  ...  ...  6j^_3^1M)  = ■ ■ IT  [B(x.)e  ^]  . 

M i=l 


k = 1 ...  M-1  O<0,  < L X.-  L 0. 

- ‘'-j-i  ^ j-i  J 


(2.257) 


Employing  the  inequalities  that  the  6 j ' s and  x ^ ' s must  satisfy, 
we  integrate  out  (x^  . . . x^^_^)  and  find 


M _Q  _ 

Xi=0i  X2-%2 


^M-l""^M-l 


V [B(x. )e  dx„  . ...  dx,  , 


(2.258) 


m ■* 


where 


(i=l 


k-1 

8.  - L x.,0 
" i=l  " 


By  considering  a series  of  linear  transformations  of 

(0,  ...  0,.  1 ) we  can  derive  the  statistics  of  the  sura  of  rela- 
1 M- 1 

tive  arrival  times  conditioned  on  M and  hence  express  the  con- 
ditional mean  as  a single  integral.  We  introduce  the  following 
two  transformations: 


(2.259) 


(2.260) 


tJm  ^ is  the  variable  that  corresponds  to  the  sum  of  the  relative 
arrival  times.  Employing  (2.259)  and  (2.260),  we  obtain  the  fol- 
lowing composite  transformations  from  (0^  ...  6j^_^)  to  ‘ ’^M-l^ 

and  the  corresponding  inverse  transformation. 
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— - - 

/”l  ^ 

Vm-i/ 

yn-J 

The  absolute  value  of  the  determinant  of  the  Jacobian  of  the 
inverse  transformation  specified  by  (2.262)  is  one.  Hence,  we 

1 

compute  p(tJj^  ...  ^ 1M)  by  making  the  substitutions  in  I 

p(6^  ...  ^ given  by  (2.258)  for  the  0 j ' s in  terms  of  the 

9i  - 

02  = t]2  ■ 2t/j^  (2.263) 

^ 9j  - T,.  - 2„._1  + ^..2 


The  constraints  on  the  tj.'s  are  obtained  by  guaranteeing  0.  > 0. 

•3  J 


j 
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> 0 


Vo  ^ 2t], 


(2,264) 


"J  2 - ^j.2 


j = 3 . . . M-1 


Hence,  the  desired  conditional  mean  may  be  expressed  as 


I?  - J J •••I 


••  I 


•••  ••• 


(2.265) 


M/N/1  Bias  Calculation 

Now  we  calculate  the  first  and  second  order  terms  in  p for 
the  asymptotic  bias  b_  „ of  the  time-contraction  algorithm  in 
the  case  of  an  M/M/1  queue.  Similar  to  our  analysis  of  the  cus- 
tomer-removal algorithm,  we  need  to  calculate  and  Oq,  the 

coefficients  of  p and  p in  an  expansion  of  z defined  by  (2.255). 

.M-1  ( ” R i 

From  the  — term  in  p(0.  ...  0„  . |M) , we  see  that  E<  E (t.  lM)f„, 
^M  ^ ( j=2  ) 

will  be  representable  as  a power  series  in  p with  terras  of  at 

least  (M-l)'st  order.  Hence,  to  compute  a.  we  need  to  calculate 

” R ^ i+1 

E L (T?\M)f|^  for  M = 2 ...  j+2,  and  collect  coefficients  of  p-^ 

i=2 


*:  *7:  . ? «> 


- , i 
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We  start  by  calculating  T2U.  From  (2,258),  p(0^  (M=2)  is 
given  by 


-Xx- 


p(9j^lM=2)  = j B(Xj^)e  . (2.266) 


V®1 


R 


Since  T2  ~ 6j^>  we  can  express  the  desired  conditional  mean 
T.J  \2  as  a double  integral  using  (2.266).  Changing  the  order  of 
integration  yields 


To\2 


2 ^ x^B(x^)e  dx^ 

x^=n 


(2.267) 


Using  the  B(x)  and  B*(S)  for  an  M/M/1  queue  specified  by  (2.196), 
we  find 


7^(2  ^ i i e Ji. 


(2.268) 


/ 


Next  we  evaluate  E 
given  by  (2.258)  as 


M 


B, 


L TjlMj  for  M=3.  p (0^^ , 62  lM=3)  is 


. / V J 


Employing  B(x)  and  B*(S)  for  an  M/M/1  queue,  we  obtain 


E f S t^1M=3|  = 7 - ^ — 

\i=2  ^ (1+p) 


6 f. 


(2.272) 


Using  Eqs.  (2.268)  and  (2.272),  we  can  calculate  the  coef- 

2 

ficients  of  the  p and  p terms  in  the  asymptotic  bias  b . a„ 

X • w • L/ 

is  derived  by  taking  the  coefficient  of  p in  an  expansion  of 

— T>  2 

(T2\2)f2‘  obtained  by  collecting  powers  of  p in 

(T2l2)f2  and  (t^  + Hence,  we  find 


1 1 
2 p 


Employing  Eq.  (2.137),  we  specify  b as 

1 • • 


(2.273) 


11  .5  12..,  3. 

2 pP  +0(p) 


(2.274) 


Hence,  the  asymptotic  bias  of  the  time-contraction  algorithm  for 
an  M/M/1  queue  contains  all  powers  of  p.  This  contrasts  with  the 
case  of  an  M/D/1  queue,  where  we  were  able  to  prove  asymptotic 
unbiasedness  of  the  algorithm. 


We  now  conclude  our  analysis  of  the  time-contraction  algor- 
ithm by  deriving  a Cramer-Rao  bound  on  the  variance  of  the  esti- 
mator for  an  M/D/1  queue.  The  form  of  the  time  contraction 
estimator  defined  in  (2.243)  and  (2.244)  tells  us  that  the  obser- 
vations which  the  algorithm  uses,  in  the  case  of  an  M/d/1  queue, 
consist  of  the  number  of  customers  served  in  each  busy  period. 
Since  the  number  of  customers  served  is  independent  from  one  busy 
period  to  another  and  for  an  M/D/1  queue  the  probability  of  M 
customers  in  a busy  period  is  given  by  Eq.  (2.110),  p(M^^  ...M^^\p) 
is  as  follows: 

M -1 

N (M.p)  -M  p 

p(M,  ,..  M^Ip)  = 7T  e . (2.275) 


Letting  y denote  the  delay  gradient,  there  is  a 1-1  corres- 
pondence between  y and  p specified  by  Eq . (2.112).  Employing 
(2.113)  ar d (2.112),  as  we  did  before  for  the  customer-addition 
and  cus toner-removal  algorithms,  we  can  derive  a Cramer-Rao  bound 
for  the  variance  of  any  unbiased  estimator  of  y that  employs 
observations  (M^^  . . . Mj^)  . 


The  bound  (2.276)  is  idenLiciil  to  that  for  the  customer 


removal  algoritltm.  This  is  to  bo  expected  since  the  maximum- 
likelihood  estimate  for  p that  follows  from  (2.275)  is  the  same 
as  that  which  follows  from  P(Y'  jp)  defined  in  (2.223),  where  Y' 
denotes  the  set  of  observations  which  the  customer-removal 
algorithm  uses  in  the  case  of  an  M/d/1  queue. 

2.11  Computational  Complexity  and  Storage  Requirement  Analysis 
for  the  Three  Algorithms 

Now  we  proceed  to  analyze,  compare,  and  contrast  the  compu- 
tational complexity  and  storage  requirements  of  the  three 
estimation  algorithms  presented.  The  procedure  for  up-dating  the 
delay-gradient  estimate  in  all  three  algorithms  assumes  the 
following  form; 


L 


^(k+i)  = ^•'ik)  ^ ^ • 


(2.277) 


is  the  delay  gradient  estimate  based  on  an  observation  inter- 


val containing  k busy  periods,  ^ is  a factor  that  renormalizes 


to  correspond  to  a component  of  the  (k+1)  busy  period 


estimate.  A denotes  the  contribution  to  the  (k+1)  busy  period 
estimate  that  comes  from  the  (k+l)-st  busy  period.  Hence,  to 
compare  the  three  algoritlims,  we  examine  the  storage  and  computa- 
tional requirements  of  forming  A . 
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In  the  customer-addition  algorithm,  A is  the  expected  addi- 
tional delay  suffered  by  the  customers  of  the  (k+l)-st  busy 
period  due  to  an  arrival  in  the  entire  (k+1)  busy  period  record. 
According  to  the  notation  defined  in  Fig.  2.1,  we  need 
(I  ...  I.  ) and  (T  ...  ...  T.  .,)  to  compute  A.  I denotes 
the  first  idle  period  at  which  an  extra  arrival  with  a service 
requirement  x can  influence  the  (ij.+l)-st  or  current  busy  period. 
Hence,  x satisfies 

X > E I,  . (2.278) 

k=nj.+l 

We  need  two  buffers  with  1^  - n^+l  storage  locations  for  the  idle 
period  and  busy  period  durations,  respectively.  We  can  general- 
ize this  argument  to  say  that  the  probability  of  x being  greater 
than  the  sum  of  idle  period  durations  is  equal  to  the  probabil- 
ity that  we  need  two  buffers  with  Nj^+1  locations.  Hence,  for  the 
case  of  an  M/G/l  queueing  system,  the  probability  of  x being 
greater  than  the  sum  of  independent  exponential  variates  with 
parameter  X is  equal  to  the  probability  of  a buffer  overflow  given 
we  have  only  storage  locations  per  buffer.  In  designing  our 
system,  we  want  to  choose  sufficiently  large  to  make  the 
probability  of  an  overflow  less  than  some  acceptable  probability 
c.  In  Eqs.  (2.44)  through  (2.46),  we  derive  the  statistics  of 
the  sum  of  k exponential  random  variables  with  parameter  X. 
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Employing  (2.46)  and  letting  p ^ Xx  we  obtain  the  following 
expression  for  the  probability  of  an  overflow  given  our  two  buf- 
fers have  N,  locations: 
b 

b ~ J 

Pr  [OverflowlN.  Storage  Location jBuffer]  =1-  e ^ L ^ . 

t>  j=o  -J' 

(2.279) 


By  applying  Taylor's  remainder  theorem  to  an  expansion  of  e^  we 
obtain 


e^  < S 
j=0 


(2.280) 


By  employing  (2.280)  we  can  upper  bound  the  probability  in  (2.279). 


Pr 


[Overflow  iNj^  Location  \ Buffer] 


(2.281) 


Since  p = Xx  and  the  range  of  interest  of  X is  0 < X < -,  we 
obtain  the  final  bound 


- ^1 
C X 1 X ^ 

Pr  [OverflowlN,  Location j Buffer}  < fpf — 


(2.282) 


We  want  to  make  the  right  side  of  the  inequality  in  (2.282)  less 

than  some  tolerable  probability  of  overflow  c.  As  an  example, 

for  X = x and  c = lO'^,  > 13  satisfies  (2.282).  Having 
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r 


1 


arrived  at  a suitable  value  of  N|^,  the  storage  requirement  of  the 
customer-addition  algorithm  is  approximately  2Nj^  + 12. 

We  now  consider  the  computational  requirements  of  forming  A . 

For  the  customer-addition  algorithm,  A is  composed  of  a contri-  | 

bution  which  represents  the  additional  system  time  resulting  from 

j 

an  arrival  in  the  (k+l)-st  busy  period  plus  the  effect  on  the  ] 

(k+l)-st  busy  period  of  an  arrival  at  any  earlier  time  in  the 
(k+1)  busy  period  record.  We  consider  first  the  contribution  due  i 

to  the  effect  of  an  extra  arrival  in  the  (k+l)-st  busy  period.  I 

i 

This  term  assumes  two  forms  depending  on  whether  all  the  service  ^ 

requirements  are  Identical.  We  employ  the  notation  defined  in  | 

Fig.  2.1  for  expressing  the  contributions  tp  D.  If  = x,  I 

the  terra  is  given  by  3 


^ (|  mV)  , (2.283) 

which  requires  four  multiplications  and  one  division  to  evaluate. 
For  the  more  general  case  qhen  x,  the  desired  contribution 

is  given  by  Eq . (2.14)  normalized  by  -t' . With  some  manipulation, 
the  component  of  A due  to  an  extra  arrival  in  the  (k+l)-st  busy 
period  is  given  by 


X 

n 


<Vi  * 


(t 


k+1 
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+ X 


M 

E 

i=2 


(2.284) 


The  above  takes  M+1  multiplications,  M+1  divisions,  and  6M-4 

additions  to  calculate.  The  T.  denote  arrival  times  relative  to 

1 

the  start  of  the  busy  period.  The  expression  in  (2.284)  is  con- 
ceptualized as  being  evaluated  as  the  busy  period  progresses. 
Hence,  we  needn't  store  all  the  service  requirements  and  arrival 
times . 

To  enumerate  the  remaining  calculations  in  forming  A we 
count  the  operations  involved  in  computing  the  effect  of  an  extra 
arrival  in  the  preceding  i^  - ^^j+l  idle  periods  and  i^^  - n^.  busy 
periods  on  the  customers  in  the  (ij^+l)-st  busy  period.  If  we 
let  i,  - n^+1  ^ n, , then  n,  corresponds  to  the  current  number 
of  idle  period  and  busy  period  durations  that  we  store.  In  terms 
of  n^,  the  remaining  operations  to  compute  A are  counted  as 


7n,  - 6 additions 

D 

6n,  - 2 multiplications 

° (2.285) 

2n,  - 1 divisions 
b 

n^^  comparisons 


We  can  upper  bound  the  number  of  operations  listed  in  (2.285) 
by  letting  n^^  = N^,  where  is  our  buffer  size  chosen  to  guar- 
antee the  probability  of  overflow  being  less  than  some  threshold 
C.  Using  our  preceding  example,  we  can  let  n^^  = 


13  for  an  c = 10 


For  the  general  case  when  the  service  requirements  are  not 


all  identical,  the  number  of  operations  to  evaluate  (2.284) 
varies  with  M.  Hence  it  is  of  interest  to  bound  the  probability 
of  M exceeding  some  value  t.  For  an  M/G/1  queue,  the  simplest 
Chebyshev  bound  yields 

Pr  (M  > t]  < ^ = . (2.286) 


Hence,  as  p -*  1 we  need  to  make  t very  large  for  the  bound  to 
give  us  any  information.  Therefore,  the  number  of  operations 
required  in  (2.284)  becomes  unbounded  as  p -*  1. 


For  the  customer-removal  algorithm  we  identify  S defined  in 
Fig.  2.5  as  the  corresponding  A.  The  number  of  storage  loca- 
tions required  to  calculate  A is  approximately  M + 8.  The  M comes 
from  the  array  of  which  we  store  as  j is  varied.  The 

number  of  operations  required  is  determined  with  reference  to  the 
flow  chart  in  Fig.  2.6. 


Minimizations  or  Comparisons 

(2.287) 

4m-3  additions 
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In  the  L ime-con  t.  rac  t ion  alj^orithm  S corr(>s()onds  to  A cTS 
before.  The  storage  requirement  to  compute  A is  8.  The  number 
of  operations  obtained  by  examination  of  Fig.  2.7  is  3M-3 
additions . 

From  the  discussion  above  we  see  that  of  all  three  algor- 
ithms, the  time-contraction  algorithm  has  the  smallest  storage 
requirement.  Unlike  the  customer-addition  and  customer-removal 
algorithms,  the  time-contraction  procedure  requires  no  buffer 
with  randomly  varying  size.  In  addition,  the  computational  load 
for  the  time-contraction  algorithm  is  the  least  of  the  three. 
Hence,  the  time-contraction  procedure  is  the  least  costly  to 
implement . 


SECTION  3 


SIMULATION  RESULTS 

In  Section  2 we  derived  the  asymptotic  bias  of  the  three 
estimation  algorithms.  Since  we  are  unable  to  calculate  the 
variance  of  the  estimators  as  a function  of  N,  the  question  of 
consistency,  whether  the  estimates  converge  asymptotically  to  the 
delay  gradient,  remains  unanswered.  In  addition,  for  an  M/D/l 
queue  we  would  like  to  know  if  our  algorithms  are  asymptotically 
efficient,  whether  they  achieve  the  Cramer-Rao  bounds  derived  in 
(2.116),  (2.224),  and  (2.276).  We  would  also  like  to  investigate 
the  robustness  of  the  customer-removal  and  time-contraction 
estimation  schemes  by  seeing  how  they  perform  for  a variety  of 
queueing  systems.  We  attempt  to  provide  answers  to  these  ques- 
tions in  the  present  section  by  presenting  the  results  of  simu- 
lating all  three  algorithms  for  an  M/D/1  queue  and  simulating  the 
time-contraction  and  customer-removal  algorithm  for  M/M/l,  D/M/1, 
and  U/M/j  queues. 

We  simulate  a single-server  queue  by  the  following  recursion 
for  successive  waiting  times: 

w . , = max  [o!+x  - 9,  01  with  ox=0.  (3.1) 

n-ri  n n n i 
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X and  0 are  random  variables  corresponding  to  the  m-th  service 
n n 

requirement  and  the  inter-arrival  time  between  the  n-th  and 
(n+l)-st  customer,  respectively.  When  oo  goes  to  zero,  this 

a 

signals  the  start  of  a new  busy  period. 


We  now  describe  the  calculations  necessary  to  evaluate  the 

statistics  of  a given  estimator.  To  derive  the  mean  and  variance 

of  a k-busy  period  estimator  we  generate  a certain  sample 

size  N of  k-busy  period  records,  processing  each  to  form  an 
o 

estimate  O',,  . . . We  compute  estimates  of  the  bias  b . and  var- 

(.k)  ,1  (K.J 

iance  cr /,  x associated  with  Dh  . as  follows: 

(k)  (k) 


(k) 


1 

S 1=1 


d: 


(k)  ,i 


bX 


(3.2) 


"2 

‘'(k) 


N, 


Ng-l 


N, 


N, 


1 ^ 2 / 1 ^ 
— ■e  n ' ^ - ~ y 

\"s  l“l  ”s  1^1 


d: 


(k)  ,i 


(3.3) 


The  measure  of  performance  that  we  use  most  often  is  the  frac- 
tional rms  error,  which  we  approximate  by  employing  our  estimates 
for  the  bias  and  variance  in  (3.2)  and  (3.3). 


(3.4) 


For  all  the  queueing  systems  under  examination,  we  present  curves 

A 

of  the  fractional  rms  error  associated  with  for  k = 10, 

100,  1000  busy  periods  with  N„  = 400  Monte  Carlo  rms  for  each 


estimate.  We  experimented  with  N , trying  N = 10,  50,  100,  200, 

o ^ 

40^^,  and  found  that  Ng  = 400  insured  from  one  to  two  significant 
figures  in  our  value  for  the  fractional  error  defined  in  (3.4). 


We  do  not  present  separate  curves  for  our  estimate  of  the 

A 

bias  b (k)  defined  in  (3.2),  since  in  many  cases  our  value  for 

A 

b(k)  was  of  comparable  size  to  the  statistical  fluctuations  in  the 
quantity.  To  make  this  notion  clearer,  (3.2)  can  be  rewritten  as 


bX/ 


+ 6 . 


(3.5) 


6 is  a zero  mean  random  variable  with 


E6^  ^ 


(3.6) 


Hence,  when  our  value  for  b(k)  is  of  the  same  size  as 


we 


can  not  possibly  measure  the  actual  bias  with  any  accuracy. 

Careful  examination  of  the  results  of  the  simulations  support  the 
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S » 


conclusion  that  in  the  cases  wliere  we  could  not  ri;jasure  h(kl 


i 


'i 


I 


accurately,  it  contributed  negli^jibly  to  the  fractional  rrns  error 
in  (3.4)  and  when  b(k)  could  be  measured  accurately,  its  con- 
tribution in  (3.4)  was  non-negligible . 

Having  defined  the  measures  we  will  use  to  compare  the 
algorithms,  we  first  discuss  the  simulation  results  for  an  M/D/1 
queue.  From  earlier  analytical  results,  we  know  the  customer- 
addition  and  time-contraction  algorithms  are  asymptotically 
unbiased  for  M/D/1  and  that  the  customer-removal  algorithm  is 
unbiased  at  least  up  to  the  third  power  of  p.  Curves  of  the  frac- 
tional rms  error  for  N = 10,  100,  1000  are  presented  for  all  three 
algorithms  in  graphs  3.1,  3.2,  and  3.3.  The  consistency  of  the 
three  estimation  procedures  is  suggested  by  the  improved  perfor- 
mance with  increasing  N.  In  graphs  3.4  and  3.5,  for  N = 100  and 
N = 1000,  respectively,  we  present  the  lower  bounds  on  fractional 
rms  error  that  follow  from  our  Cramer-Rao  bounds  for  the  variance 
of  the  estimators,  together  with  the  simulation  results.  The 
closeness  of  the  simulation  curves  to  their  respective  lower 
bounds  suggest  that  all  three  algorithms  are  asymptotically 
efficient  for  an  M/D/1  queue.  Graphs  3.4  and  3.5  show  that  as  a 
function  of  p the  fractional  rms  error  for  the  customer-addition 
algorithm  remains  below  that  of  the  time-contraction  and  customer- 
removal  algorithms,  and  hence  we  conclude  that  it  performs  the 
best  of  the  three  procedures  for  an  M/d/1  queue.  The  customer- 
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removal  and  time-contraction  algorithm's  fractional  rins  error 

curves  are  quite  close,  with  the  customer-removal  algorithm's 

performance  being  slightly  worse.  The  bias  results  for  the  three 

algorithms  are  displayed  in  Tables  3-1  and  3-2.  We  list  our 

estimate  for  the  relative  bias  b(k)/— --  together  with  the  size  of 

^ a(k)/,,'Ng 

statistical  fluctuations  in  the  relative  bias,  . 

0 u 

Seemingly  anomalous  behavior  of  the  relative  bias  is  often 
explainable  by  looking  at  the  size  of  the  statistical  fluctuations 
in  the  quantity. 

We  now  proceed  to  demonstrate  that  the  behavior  of  the 
algorithms  for  p near  zero  indicated  by  graphs  3.1,  3.2,  and  3.3 
is  reasonable.  We  examine  all  three  algorithms  for  p — 0 by 
keeping  x constant  and  considering  X — 0.  As  p is  made  arbitrarily 
close  to  zero,  the  density  for  M^,  the  number  of  customers  served 
in  the  i-th  busy  period,  becomes  an  impulse  at  1.  Hence,  as  p -»  0, 
the  queueing  record  approaches  a set  of  single  customer  busy  per- 
iods. In  addition,  since  I , the  average  size  of  the  k'th  idle 

K 

period  is  the  idle  periods  become  unbounded  as  p -*  0.  Hence, 
with  high  probability  the  service  requirement  of  the  additional 
customer  x will  be  « Applying  (2.14),  (2.15),  and  (2.17) 

through  (2.19),  we  find  that  with  high  probability  E(6S\tcTj^, 
Queueing  Record)  and  E(6S\tcIj^,  Queueing  Record)  will  be  given  as 


E(6SltfTj^,  C j jueing  Record)  = x + y x 


(3.7) 


■*r  ’ir'  .f  '■  — — ' 


ft  ♦ . , 


TABLE  3-2 


RELATIVE  BIAS  FOR  CUSTOMER-REMOVAL  AND 
TIME-CONTRACTION  ALGORITHM  (M/D/1  QUEUE) 


p 

N 

‘-C.R. 

CJAi 

CN).C.R. 

\.C. 

a^,  /vN 

(N),T.C. 

6D 

^.X 

cx 

10 

1.6  X 10"^ 

8.7  X 10'^ 

8.8  X 10"^ 

6.8  X 10"^ 

.1 

100 

1.6  X 10"^ 

2.2  X lO"^ 

7.8  X 10"^ 

1.9  X lO"^ 

1000 

4.7  X 10“^ 

7.9  X 10“^ 

2.8  X 10"^ 

7.3  X 10~^ 

10 

1.2  X lO"^ 

1.7  X 10‘^ 

-2 

5.4  X 10 

1.4  X 10“^ 

.3 

100 

6.3  X lO"^ 

5.1  X 10’^ 

-3.2  X 10’^ 

4.7  X 10“^ 

1000 

5.0  X lO"^ 

1.9  X 10'^ 

2.8  X 10’^ 

1.8  X lO'^ 

10 

3.1  X 10"^ 

2.4  X 10"^ 

-3.0  X 10"^ 

2.3  X 10’^ 

.5 

100 

-1.4  X 10"^ 

9.8  X 10“^ 

-1.7  X 10’^ 

9.2  X lO"^ 

1000 

-9.0  X 10"^ 

3.5  X lO"^ 

-1.4  X lO'^ 

3.3  X lO'^ 

I 

10 

-1.7  X 10'^ 

2.8  X 10~^ 

-1 

-1.8  X 10 

2.8  X 10'^ 

100 

-5.3  X 10"^ 

1.6  X 10'^ 

-5.7  X 10"^ 

1.6  X lO"^ 

1 

1000 

3.2  X lO"'^ 

6.3  X lO'^ 

-4.6  X 10  ^ 

6.2  X lO"^ 
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j,  . 7- 


V f 


r 


1 X 


E(6S\tcI.  , Queueing  Record)  = x + ^ — 


(3.8) 


N-1 

Using  the  fact  that  T„  ~ Nx  + L L for  p -*  0 and  substituting 

^ k=l 

(3.7)  and  (3.8)  into  (2.5)  and  (2.6)  to  form  the  delay  gradient 

A 

estimator,  we  obtain  the  following  expression  for  : 


''(N)  ~ " 


(2  - — ) — x^ 

1 N-1 

^ N ^ 

N k=l 


(3.9) 


If  we  let  p approach  zero,  we  can  make  the  probability  that 

1 N-1 

— S I.  is  greater  than  any  given  number  converging  to  1 (since 

N k=l  ^ 

I.  = r)  . Hence,  in  the  limit  as  p -«  0 the  second  term  in  (3.9) 

K A 

yields  a zero  contribution  and  ~ Therefore,  not  only  is 

Ai 

the  customer-addition  algorithm  unbiased  as  p -»  0,  but  Var  — 0 

for  p 0. 

Since  as  seen  in  Eqs.  (2.123)  and  (2.242),  the  customer- 
removal  and  time-contraction  algorithms  have  a similar  structure, 
we  can  examine  the  behavior  of  both  concurrently  for  p -»  0.  As 
p 0,  ~ 1 and  each  customer  has  a zero  waiting  time.  and 

in  Eqs.  (2.123)  and  (2.242)  become  zero,  and  both  estimators 
assume  the  following  form: 


159 


■ V 


L X, 


(3.10) 


X,  ' denotes  the  service  requirement  of  the  first  customer  in  the 

j-th  busy  period.  Hence,  for  p - 0 ~ x and  Var  , 

2 

where  a,  denotes  the  variance  of  the  service  time  density.  Since 
b 


from  Eq.  (2.24)  it  follows  that 


= X for  a general  single- 


server queue,  the  customer-addition  and  time-contraction  algor- 
ithms will  be  unbiased  near  p = 0,  but  the  estimators  will  have  a 
2 

variance  cTj^/N.  Hence,  since  = 0 for  an  M/D/1  queue,  we  expect 
graphs  3.2  and  3.3  of  the  fractional  rms  error  to  pass  through 
the  origin. 

Before  examining  the  robustness  of  the  time-contraction  and 
customer-removal  algorithms,  we  review  the  theoretical  relations 
for  the  delay  gradient  in  the  case  of  M/G/1  and  G/M/1  queues. 
Employing  (2.130)  for  in  Eq . (2.22),  we  obtain  the  following 
relation  for  an  M/G/1  queue: 


bD  - , . 


(1  + C^) 


(1  -p)‘ 


(3.11) 


c2  A Ik 

^b  - -2  ■ 

X 


(3.12) 


For  an  M/D/1  queue  C,  = 0,  while  for  an  M/M/1  queue  C,  = 
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In  a G/M/1  queue  we  have  a general  inter-arrival  time  density 
A(x)  with  one-sided  Laplace  transform  A*(s),  and  an  exponentially 
distributed  service  requirement  with  parameter  p.  The  waiting 
time  W is  specified  by 

_ = q (3.13) 

W p(l  - a)  ’ 

where  a solves  the  nonlinear  equation 


Relation  (3.14)  is  expressed  as 

» JL  Z. 

a = e P e^  . (3.18) 

Differentiating  (3.18)  with  respect  to  p we  find 


ba  ^ g (tng) 
bp  gtng  + (1-g) 


(3.19) 


For  a U/M/1  queue,  the  inter-arrival  times  are  uniformly  dis- 
tributed. A(x)  and  A*(s)  are  given  by 


A(x)  = j [lJ_^(x)  - U_^(x  - |-)] 

(3.20) 

/ 2 \ 

A*(s)  =1  i[l  - e^  '1 

(3.21) 

U_j(x)  denotes  the  unit  step  function.  Relation  (3.14)  may  be 


expressed  as 


g 


- e 


2(1  -g) 


(3.22) 


Differentiating  (3.22)  with  respect  to  p we  obtain  — as 

op 


bo 

bp 


1 - e 


P 1 4 2 

V D 


(1-g) 


:(l 


2g  + e ” / 


(3.23) 
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To  construct  ^ as  a function  of  p for  a G/M/1  queue,  we 
Ca 

vary  p,  solving  for  the  appropriate  a by  employing  the  fixed- 
point  iteration  method  on  Eq . (3.14) 

= A*(p  - po^)  (3.24) 

We  select  some  starting  0 < Oq  < 1 and  apply  the  above  iteration 

until  a .1  = a to  the  desired  accuracy. 
n+1  n 

Now  we  examine  the  robustness  properties  of  the  customer- 
removal  and  time-contraction  algorithms  by  comparing  their  per- 
formance for  M/M/1,  U/M/1,  and  D/M/1  queues.  The  simulation 
results  for  the  fractional  rms  error  of  the  two  algorithms  for 
M/M/1,  U/M/1,  and  D/M/1  queues  are  presented  in  graphs  3.6  through 
3.11,  respectively.  The  results  for  the  relative  bias  are  dis- 
played in  Tables  3-3  through  3-5.  The  time-contraction  algorithm 
performs  slightly  better  than  the  customer-removal  algorithm  for 
M/M/1  and  both  perform  nearly  the  same  for  a U/M/1  queue.  The 
only  dramatic  differences  occur  for  the  D/M/1  queue.  Both  the 
relative  bias  and  fractional  rms  error  show  the  time-contraction 
algorithm  performing  better  than  the  customer-removal  procedure. 
This  result  is  reasonable  since  for  a D/M/1  queue  the  time- 
contraction  algorithm  simulates  a change  in  arrival  rate  in  the 
exact  way  dictated  by  the  structure  of  the  queue. 
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TABLE  3-3 


RELATIVE  BIAS  FOR  CUSTOMER- REMOVAL  AND 
TIME-CONTRACTION  ALGORITHM  (M/M/1  QUEUE) 


N)  .C.R 


5D 


bX 

bX 

bX 

— 

10 

-1.5  X lO"^ 

1.8  X lO'^ 

-1.4  X lO"^ 

1.8  X lO"^ 

100 

1.3  X lO"^ 

7.4  X lO"^ 

2.3  X lO"^ 

6.7  X 10"^ 

1000 

7.2  X 10“^ 

2.5  X 10“^ 

8.5  X lO"^ 

2.4  X lO"^ 

10 

-2.3  X lO"^ 

2.4  X lO'^ 

-2.2  X lO"^ 

2.4  X lO'^ 

100 

-1.5  X 10"^ 

1.3  X 10“^ 

2.6  X lO"^ 

1.2  X lO"^ 

1000 

2.1  X 10'^ 

4.4  X 10"^ 

2.4  X 10'^ 

4.1  X lO"^ 

10 

-3.0  X 10'^ 

3.2  X 10"^ 

-3.0  X lO"^ 

3.2  X 10‘^ 

100 

-2.3  X 10"^ 

1.9  X 10"^ 

-2.8  X 10"^ 

1.7  X 10"^ 

1000 

3.2  X lO"^ 

6.4  X 10"^ 

3.6  X 10"^ 

6.3  X 10“^ 

10 

-4.4  X 10"^ 

3.3  X 10"^ 

-4.5  X 10'^ 

3.1  X lO"^ 

100 

-3.0  X 10'^ 

2.8  X 10"^ 

-2.5  X 10"^ 

2.7  X 10"^ 

1000 

6.4  X 10'^ 

1.1  X 10"^ 

8.8  X lO'^ 

1.1  X 10"^ 

TABLE  3-4 


RELATIVE  BIAS  FOR  CUSTOMER- REMOVAL  AND 
TIME- CONTRACT ION  ALGORITHM  (D/M/1  QUEUE) 


TABLE  3-5 

RELATIVE  BIAS  FOR  CUSTOMER- REMOVAL  AND 
TIME-CONTRACTION  ALGORITHM  (U/M/1  QUEUE) 


N)  .T.C 


e>D 


bX 

6X 

5X 

1 

1 

10 

-1.4  X lo"^ 

1.7  X 10'^ 

-1.3  X lO"^ 

1.7  X lO"^ 

100 

-1.2  X lo"^ 

6.4  X lO"^ 

1.5  X lO"^ 

6.1  X 10"^ 

1 

1000 

i 

5.2  X 10“^ 

2.1  X 10’^ 

7.1  X 10‘^ 

2.0  X 10“^ 

10 

-2.1  X 10"^ 

2.1  X lO"^ 

-1.9  X 10"^ 

2.3  X 10“^ 

100 

-3.5  X lO"^ 

1.1  X 10"^ 

-2.7  X 10"^ 

1.1  X 10'^ 

1000 

-1.3  X 10"^ 

3.6  X 10"^ 

3.1  X lO'^ 

2.8  X 10'^ 

10 

-3.5  X 10"^ 

2.7  X 10'^ 

-2  9 X 10"^ 

3.1  X 10"^ 

100 

-9.7  X 10"^ 

1.5  X 10'^ 

-1.2  X 10'^ 

1.6  X 10"^ 

1000 

-4.3  X lO"^ 

6.3  X lO"^ 

4.5  X lO"^ 

6.4  X lO"^ 

10 

-4.7  X lO"^ 

3.2  X 10'^ 

-4.1  X 10'^ 

3.5  X 10"^ 

100 

-1.4  X 10“^ 

2.4  X 10“^ 

-3.7  X 10'^ 

2.5  X 10'^ 

1000 

-4.2  X 10"^ 

9.4  X 10’^ 

5.7  X 10"^ 

9.8  X 10’^ 

In  examining  the  fractional  rms  error  curves  of  the 
customer-removal  and  time-contraction  algorithms  for  M/M/1, 

U/M/1,  and  D/M/1  queues,  we  first  note  that  the  behavior  of  each 
algorithm  is  nearly  the  same  for  all  three  queues  when  p is  near  0. 
This  follows  from  our  demonstration  that  the  estimator  correspond- 
ing to  both  algorithms  approaches  Eq.  (3.10)  as  p — 0.  Since 

for  all  our  simulations  we  take  x ~ ^ ~ 1»  the  variance  of  the 

2 

service  time  density  is  = 1.  Hence,  our  N busy  period  curve 

for  the  fractional  rms  error  for  each  algorithm  should  cross  the 

axis  at  approximately  -3.  The  simulation  results  show  this  most 

VN 

clearly  in  the  curves  for  the  M/M/1  and  U/M/1  queues.  The  points 
derived  from  the  simulation  of  the  D/M/1  queue  are  more  scattered 
than  those  for  the  other  two  queues,  making  the  interpretation  of 
behavior  near  p = 0 uncertain. 

For  each  queue  the  performance  of  the  estimation  algorithms 
degrades  steadily  with  increasing  p.  Since  the  queue  becomes 
non-stationary  at  p = 1 , it  is  reasonable  that  the  fractional 
rms  error  should  become  unbounded  at  p = 1 . We  can  summarize  the 
performance  of  the  time-contraction  and  customer-removal  pro- 
cedures for  each  queue  by  listing  least  upper  bounds  on  the  frac- 
tional rms  error  for  p in  the  interval  [0,.7] . For  the  N = 100 
and  N = 1000  busy  period  estimators,  tliese  least  upper  bounds 
are  listed  in  Table  3-6. 
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TABLE  3-6 


N = 1000 

N 100 

M/M/1 

D/M/1 

U/M/1 

M/M/1 

D/M/1 

U/M/1 

Cus  tomer-Removal 
Algorithm 

.21 

.25 

.20 

.56 

.46 

.47 

Time-Contraction 

Algorithm 

.22 

.19 

.20 

.54 

.45 

.51 

m 


We  now  pose  the  question  of  how  long  the  observation  period 

T must  be  to  contain  an  average  number  of  N busy  periods  and 
E 

hence  achieve  the  fractional  rms  errors  reported.  A record  of 
length  T contains  an  average  number  of  customers  XT  . From 
queueing  theory  we  can  derive  an  expression  for  M,  the  average 
number  of  customers  served/busy  period.  Hence,  \T^/M  will  be 
the  average  number  of  busy  periods/T^  sec.  observation  time. 
Appealing  to  M/G/1  queueing  theory,  M is  given  by  Eq . (2.93)  and 
we  can  obtain  the  following  condition  on  T so  the  average  number 
of  busy  periods  contained  in  the  observation  period  is  greater 
than  or  equal  to  N. 


T 


E 


X N 

- p(l  - p) 


(3.25) 


We  now  interpret  (3.25)  for  the  case  of  queues  that  occur 
in  computer  networks.  In  this  situation,  the  customers  are 
messages  with  a certain  number  of  bits  1.  The  service  require- 
ment is  the  time  needed  to  transmit  the  message  t/c,  where  c is 
the  capacity  of  the  communications  link  in  bits/sec.  Hence  x 
is  given  by 


X 


(3.26) 


For  representative  purposes,  values  of  I,  c,  and  x for  the  Arpa 
network  are  listed  below 
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Jx  . 


J 


t = 1000  bits 

c “ 50  X 10^  bits/sec  (3.27) 

1 

X “ sec. 

If  p is  unrestricted,  (3.25)  will  require  an  unbounded  observa- 
tion time  as  p is  made  arbitrarily  close  to  0 or  1 . Hence,  we 
hypothesize  that  p is  restricted  to  the  interval  [.1,  .9]. 
Employing  the  value  in  (3.27)  for  x and  letting  pc  [.1,  .9],  for 
the  observation  interval  T„  to  contain  an  average  of  N busy  per- 
iods  we  must  have 


T„  > .222  N.  (3.28) 

E — 

For  N = 1000  Eq.  (3.28)  yields  T^.  > 222  sec.  or  T^,  > 3.7  min. 
and  for  N ••  100,  T„  > 22.2  sec.  or  T_  > .37  minutes. 

Cl  *”  Cl 
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SECTION  4 
CONCLUSION 


Hence,  we  have  three  estimation  algorithms  that  appear  from 
simulation  and  analytical  results  to  be  asymptotically  unbiased, 
consistent,  and  asymf)toLically  efficient  in  the  case  of  an  M/d/1 
queue.  We  proved  asymptotic  unbiasedness  for  the  time-contraction 
and  customer-addition  algorithms  in  the  M/D/l  case.  We  showed 
that  the  asymptotic  bias  for  the  customer-removal  algorithm, 
expressed  as  a power  series  in  p,  only  contains  terms  of  third- 
order  or  higher  (for  an  M/D/1  queue).  However,  the  closeness  of 
the  customer-removal  simulation  results  to  those  of  the  time- 
contraction  algorithm  suggest  that  the  customer-removal  pro- 
cedure is  asymptotically  unbiased.  The  consistency  and  asymp- 
totic efficiency  of  the  three  algorithms  follows  from  comparing 
the  fractional  rms  errors  derived  from  simulating  the  procedures 
on  an  M/D/1  queue  with  lower  bounds  on  fractional  rms  errors 
derived  from  Cramer-Rao  bounds  on  the  variance  of  any  unbiased 
estimator . 

In  evaluating  the  most  promising  algorithm  as  far  as  robust- 
ness, computational  complexity,  and  storage  requirements,  we  can 
only  choose  between  the  customer-removal  and  time-contraction 
algorithms,  since  the  customer-addition  algorithm  as  formulated 
is  only  applicable  to  queues  where  all  customers  have  the  same 


.x:  . . 


service  requirement.  The  customer-removal  algorithm  requires  on 
2 2 

the  order  of  M minimization,  M additions,  and  M storage  loca- 
tions to  form  the  quantity  needed  to  up-date  the  estimate  at 
the  end  of  the  current  busy  period,  where  M is  the  number  of  cus- 
tomers served  in  the  most  recent  busy  period.  The  time-contraction 
algorithm  requires  eight  storage  locations  and  on  the  order  of  M 
additions.  In  examining  the  performance  of  the  two  algorithms 
for  M/D/1,  M/M/1,  D/M/1,  and  U/M/1  queues,  the  only  dramatic  dif- 
ference occurs  in  the  case  of  a D/M/l  queue,  where  the  customer- 
removal  algorithm  behaves  worse  as  reflected  in  the  size  of  the 
relative  bias  and  the  scatter  of  simulation  points  for  the  frac- 
tional rms  error.  Analytical  results  show  that  both  algorithms 
are  not  asymptotically  unbiased  for  an  M/M/1  queue.  Hence,  our 
• results  suggest  the  time-contraction  algorithm  as  the  best  candi- 

date since  it  appears  to  be  at  least  as  robust  as  the  customer- 
1$  removal  procedure,  while  having  considerably  less  computational 

and  storage  requirements  . 

Whatever  the  estimation  procedures  we  employ,  there  is  a 
trade-off  between  observation  time  and  the  accuracy  of  our  estimates 
that  needs  to  be  investigated  further.  If  we  up-date  our  routing 
• variables  every  T^  seconds,  it  would  be  desirable  that  our 

estimates  should  converge  in  approximately  T^  seconds.  On  the 
basis  of  this  criterion,  and  using  the  result  for  the  Arpa  net 
derived  in  Eq . (3.28),  if  we  up-date  our  routing  variables  every 
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= 7.4  minutes,  we  will  obtain  estimates  with  performance  at 
least  as  good  as  the  N = 1000  busy  period  estimator.  Equation 
(3.28)  assumes  that  p is  confined  to  [.1,  .9].  When  the  link  is 
very  free  or  very  loaded  the  estimates  will  not  be  as  accurate. 
However,  this  is  ameliorated  by  the  fact  that  in  the  loaded  case 
we  will  try  to  take  traffic  away  and  in  the  lightly  loaded  case 
we  will  try  to  add  additional  traffic.  In  the  Arpa-network , 
delay  estimates  are  exchanged  between  nodes  and  new  routes  may  be 
determined  every  ^ second.  Hence,  we  may  not  be  able  to  up-date 
the  routing  as  often  as  we  like  if  we  expect  to  do  the  estima- 
tion of  the  delay  gradients  with  reasonable  accuracy.  However, 
since  in  a quasi-static  procedure  the  changes  in  the  routing 
variables  are  likely  to  be  small,  highly  accurate  estimates  may 
be  unnecessary. 

One  possibility  for  future  work  is  in  investigating  the 
convergence  of  the  quasi-static  routing  algorithms  presented  in 
Section  1 given  that  they  use  estimates,  instead  of  exact  values, 
for  the  marginal  delays.  The  most  ambitious  approach  would  be 
to  simulate  the  entire  communication  net  with  queues  at  each  link 
and  apply  the  estimation  algorithms  and  the  different  "control" 
strategies  for  adjusting  the  routing  variables.  A second  approach 
might  be  to  approximate  the  estimates  as  some  random  variable  we 
generate  in  the  computer  with  a mean  and  variance  derived  from 
simulation  results  and  see  whether  the  quasi-static  routing  algor- 
ithms still  converge  to  the  optimal  flow  pattern. 
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APPENDIX  A 


VERIFICATION  TilAT  (2.117)  IS  THE  SOLUTION  OF 
DIFFERENCE  EQUATIONS  (2.176) 


The  proof  that  (2.177)  is  the  solution  of  the  difference 
equations  specified  by  (2.176)  follows  from  substituting  the 
solution  into  (2.176)  and  checking  for  equality.  The  following 
two  identities  are  useful  at  certain  steps  in  the  proof: 


(q^X-i)*^'  = 0 


(A.l) 


L ( ^)(-l)‘^'q'  « 0 (for  k > 1)  (A. 2) 

q'-l  ^ 

Formula  (A.l)  follows  from  using  the  binomial  representation  of 

(1-x)  and  substituting  x = 1.  Formula  (A. 2)  comes  from  similarly 
d k 

expanding  — (1-x)  and  evaluating  it  at  x = 1 for  k > 1. 

Substituting  our  expressions  for  and  given  by 

1 q 

(2.177)  into  Eq . (2.176),  cancelling  a common  factor  of  x^  ^/(q-1) 
that  appears  on  both  sides,  and  making  Lhe  change  of  variables 
j •=  n-2-i,  we  are  left  with  the  following  relation  which  we  must 
prove. 


n 


n-q-1 


(n-q)! 


n-q-1 

L 

j=0 


Cn.-2-..i) 


q(n-l-q-j) ! (j+1) ! 


(n-D- 


q(n-q) ! 


(A. 3) 
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Letting  q = n-k,  we  can  recognize  the  binomial  coefficient 
and  multiplying  both  sides  by  k! (n-k)  we  are  left  with 


k k-1  k i 

n - kn  = E (. .i)(n-2-j)(n-l)J  + 1. 

j=0  J 


(A. 4) 


Relation  (A. 4)  must  be  true  for  k = 0 ...  n-1  with  n > 1. 

By  inspection,  the  relation  is  true  for  k=0  and  k=l . We  prove  (A. 4) 
for  k > 1 by  matching  the  coefficients  of  powers  of  n on  both 
sides  of  the  equation.  We  substitute  for  (n-1)^  the  representa- 
tion given  by  the  binomial  formula  and  then  change  the  order  of 
summation,  obtaining  the  following  expression  for  the  right  side 
of  (A. 4): 
k-1 


k 

E 

4-1 


j=4-l 


n 


k-1 
+ E 
4=0 


k-1 

E 


n^  + 1. 


(A. 5) 


The  coefficient  for  n in  (A. 5)  is  given  by 
j=k-l  ^ • 


(A. 6) 


k-1 

The  coefficient  of  n is  given  by 


j=k-2 


'j+l''j-k+2' 


j=k-l 
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^j+1^  "j-k+1' 


(A. 7) 


The  coefficient  of  is  obtained  by  considering  the  second  and 
third  terms  in  (A. 5). 


V (,+T)(-5)(-l)-^'^^(2+j)  + 1 

j=0  J -• 


(A. 8) 


Making  the  change  of  variable  j =q'-l  and  then  employing  (A,l) 
and  (A. 2),  since  we  consider  k > 1,  we  obtain  the  following 
expression  from  (A. 8): 


E (q'+l)  + 1 

q‘-l 


E (^X-l)'’’  + E ( ^)(-l)'’'q'  =0, 
q’=0  ^ q'=l  ‘I 

(A. 9) 


■*  ' 


The  final  step  in  the  proof  is  to  show  that  the  coefficients  of 
n“^  for  -1  = 1 ...  k-2  are  zero.  From  (A. 5)  the  desired  coefficient 
of  n^  is  expressed  as 


(A. 10) 


Making  the  change  of  variables  j ~q'  +-t-l,  breaking  the  second 
term  into  two  summations,  and  regrouping  the  binomial  coefficient 
terms  in  one  of  the  new  summations,  we  obtain 


From  (A. 2),  for  k-t  > 1 or  1 < t < k-2 , we  recognize  the  last 
term  in  (A. 11)  as  zero.  Grouping  the  first  two  summations  in 
(A. 11)  together,  we  are  left  with 


Substituting  (A. 15)  into  (A. 14),  adding  and  subtracting  1 in  the 
summation,  we  are  left  with 


msBSSssaia 


(V)(-i)'^'  - 1 

^ ^ (q'=0 


(A. 16) 


But  Eq . (A. 16)  is  zero,  since  by  (A.l)  the  summation  inside  the 
brackets  is  zero. 


184 


r 


i 


APPENDIX  B 

CALCULATION  OF  C^W,  C^  \4  FOR  AN  M/D/1  QUEUE 


— 1 —1 

We  outline  here  the  calculation  of  C3\4, 

and  C^l4. 

C3  IS  given  by 

C3  = min  ix.zj  , 

(B.l) 

where 

z = min  [0)2.003]  . 

(B.2) 

The  joint  density  of  ^2  and  conditioned  on  M-3  is  as  follows; 


p(co_,co^  1M=3)  = 2 

^ ^ 3x 


0 < u;i2  < X 

0 < W3  < X + 


(B.3) 


Hence,  the  distribution  function  of  z is 


= 1 - J J 


z \M=3 

0 < T < X 


2 

3x‘ 


X + 


C02 


duJod 


u)2=t  u)3=t 


Differentiating  (B.4)  with  respect  to  t we  obtain 


fz|M-3<">  (2X-0 


0 < T < X 
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(B.4) 


(B.5) 


j 


«v 


^ I ^ • 


Since  (2.168)  C^\3  is  computed  as 


(B.6) 


We  next  consider  the  calculation  of  0^14.  and  the 

appropriate  z are  defined  in  (B.i)  and  (B.2).  The  joint  density 
of  (0^2, conditioned  on  M=4  is  specified  by  (2.169)  as 


X + w„ 


p(o;i2,uJ2  1M=4)  = — j J du)^  . 


(B.7) 


0 < 0)2  < X 

0 <00^  < X + CO2 


^4  = 0 


The  distribution  function  of  z is  found  as 


1 - 77  J I (X+1C3)  dcojdc^ 


z|M-4'" 
0 < T < X 


0)2*7  01^*7 


(B.8) 


Differentiating  (B.8)  with  respect  to  7,  ^2lM=4^^^  given  by 


0 < 7 < X 


(B.9) 


Noting  again  that  ^ (2.168)  0^14  is  calculated  as 


■ J"  ’■''zlM-4l’’>  ■ 2 

7=0 


(B.IO) 


\ ft  j>~..  • , 


Our  final  calculation  is  of  c]^14.  is  defined  by 


min  {x,z1  , 


where 


z = min  [c02,W3,w^]  • 


(B.ll) 

(B.12) 


From  the  joint  density  of  (0^,003,00^)  conditioned  on  M-4  given  in 
(2.190),  we  express  the  distribution  function  of  z as 


F 

0 


z\M=4^'^^ 

< T < X 


1 


X + X + CO3 

I I dto^du;3du)2  • 

oJ3=t 


(B.13) 


Differentiating  (B.13)  with  respect  to  t , we  obtain 


^ - 3xr  + I r^l  . 

0 < T < X 

Since  F^,j^^^(x)  = 1,  by  (2.168)  C ^ W is  computed  as 


(B.14) 
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